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ABSTRACT 


The  objective  of  this  research  is  to  develop  a  constitutive  theory  for  granular 
material  from  a  micromechanics  point  of  view  so  that  the  effects  of  micro-structure  are 
considered.  This  type  of  micro-structural  based  constitutive  theory  is  useful  in  many 
fields  of  studies  such  as  soii  mechanics,  powder  mechanics,  composite  mechanics 
and  ceramic  mechanics.  The  specific  efforts  are  focused  in  the  foiiowing  four 
different  areas:  1)  descriptions  of  micro-structure,  2)  micro-macro  relationship,  3) 
classes  of  micromechanics  constitutive  theory,  and  4)  contact  law  of  the  inter-particle 
binder.  These  four  areas  are  the  fundamental  elements  to  the  construction  of  a 
micromechanics  theory  for  granular  media. 

Previous  studies  have  shown  that  at  large  deformation,  the  strain  exhibits  non¬ 
uniformity  within  the  material  sample  which  contradicts  to  the  uniform  strain 
assumption  used  in  the  past  theories.  The  non-uniformity  of  strain  has  been  a  major 
obstacle  for  constitutive  modelling.  Therefore,  particular  attention  will  be  given  to  the 
description  of  heterogeneity  in  microstructure  and  the  methodology  of  integrating  the 
microstructure  description  in  relating  the  micro-macro  mechanical  behavior.  The 
micromechanjcs  theory  will  lead  to  a  model  that  can  capture  the  features  due  to 
effects  of  microstructure  and  discrete  nature  of  granular  particles. 

A  related  aspect  in  the  development  of  such  a  micromechanics  model  is  to 
find  a  link  between  discrete  system  and  continuum  system.  The  results  of  this 
research  will  provide  basic  understanding  in  different  classes  of  continue  suitable  for 
granular  materials  and  in  the  limitations  of  the  classic  continuum  in  the  usual 
continuum  theory.  The  developed  theory  will  be  evaluated  by  special  features  that 
exist  only  in  granular  material,  such  as  the  dispersion  and  decay  of  stress  wave 
propagating  in  granular  media.  The  nature  of  this  investigation  is  focused  on 
theoretical  development.  The  developed  theory  is  evaluated  by  experimental  and 
computer  simulation  results. 
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CHAPTERl 

Summary  of  The  Project 


1.1  Research  Objectives 

The  general  objective  of  this  research  is  to  develop  a  constitutive  model  for  a  packing 
of  granules.  The  effort  is  focused  in  four  different  areas:  1)  statistical  descriptions  of  micro¬ 
structure,  2)  micro-macro  relationship,  3)  classes  of  micromechanics  constitutive  theory,  and 
4)  contact  law  of  the  inter-particle  binder.  These  four  areas  are  essential  to  the  construction  of 
a  micromechanics  theory  for  particulate  media. 

1.2  Accomplishments 

Accomplishments  of  the  present  project  is  reviewed  in  four  different  areas;  1) 
descriptions  of  micro-structure,  2)  micro-macro  relationship,  3)  classes  of  micromechanics 
constitutive  theory,  and  4)  contact  law  of  the  inter-particle  binder. 

1.2.1.  Descriptions  of  Micro-structure 

A  proper  description  of  microstructure  of  granular  material  is  imperative  to  the 
modelling  of  its  micro-scale  behavior.  Three  approaches  of  microstructure  description  have 
been  used.  They  are: 

1)  particle  pair  (branch)  model, 

2)  particle  group  (cell  or  micro-element)  model,  and 

3)  models  with  two-point  probability  function. 

In  the  particle  pair  model,  the  configuration  of  each  particle  pair  is  represented  by  a 
branch.  The  microstructure  of  granular  material  is  represented  by  a  distribution  function  of  the 
lengths  and  orientations  of  all  branches.  The  orientational  distribution  is  usually  expressed  in 
terms  of  a  fabric  tensor  which  represents  the  anisotropic  properties  of  a  granular  media. 

Along  this  line,  we  studied  the  mechanical  properties  in  relation  to  the  material  symmetry 
represented  by  the  fabric  tensor  (Ref  13,  9,  12,  14,  and  15  in  section  1.3.3). 

A  more  detailed  description  for  microstructure  is  to  divide  the  material  into  cells  (e.g. 
Voronoi  tessellations  or  Delaunay  tessellations).  Each  cell  represents  a  group  of  particles. 
Configuration  of  each  particle  group  is  described  by  a  local  fabric  tensor.  The  microstructure 
of  granular  material  is  represented  by  a  distribution  function  of  the  local  fabric  tensors.  We 
utilize  this  more  complex  description  to  account  for  the  variation  of  microstructure^or  different 
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cells  (i.e.,  the  heterogeneity  of  material).  It  was  found  that  the  effects  of  heterogeneity  are 
significant  on  plastic  deformation  (Ref.  3,  4,  5,  7  in  section  1.3.3). 

To  further  describe  the  nature  of  heterogeneity,  we  use  the  two-point  probability 
function.  A  new  formulation  is  derived  that  determines  the  two-point  probability  function 
directly  from  the  statistical  measures  of  stereological  chord  length  and  porosity.  Two-point 
probability  function  is  a  useful  and  practical  descriptor  for  the  microstructure  of  a  general 
random  multi-phase  material.  Methods  are  developed  (24,  30)  for  predicting  mechanical 
property  based  on  the  two-point  probability  function.  Compared  to  the  conventional  effective 
medium  theory,  the  two-point  probability  method  provides  a  viable  approach  to  account  for  the 
interactions  among  phases  of  material.  Work  continues  in  the  area  of  characterizing 
heterogeneity  by  combining  the  concepts  of  (2)  and  (3).  The  objective  is  to  search  for  a 
simple  approach  that  captures  the  essential  information  of  microstructure. 

1.2.2.  Micro-macro  Relationship 

Four  methods  have  been  used  to  relate  micro-macro  quantities. 

1)  kinematic  hypothesis  :  the  particle  displacements  are  determined  from  macro-scale 
strain. 

2)  static  hypothesis  :  the  contact  forces  are  determined  from  macro-scale  stress. 

3)  self-consistent :  the  particle  displacements  are  determined  from  the  micro-element 
strain;  the  micro-element  strain  is  determined  from  macro-scale  strain  based  on 
Eshelby  method. 

4)  statistical  method  :  the  particle  displacements  are  determined  from  the  micro-element 
strain;  the  micro-element  strain  is  determined  from  macro-scale  strain  based  on  a 
statistical  distribution  function. 

In  the  kinematic  hypothesis,  the  uniform  strain  assumption  is  implied.  Therefore,  the 
derived  constitutive  constants  are  applicable  only  to  small  strain  conditions  where  the  strain  is 
relatively  uniform  in  a  sample.  At  the  level  of  large  deformation,  the  strain  becomes  highly 
non-uniform  due  to  particle  sliding  and  particle  rotation.  When  the  strain  field  is 
heterogeneous,  it  is  found  that  the  overall  behavior  is  no  longer  a  straight  volume  average  of 
the  local  behavior  (Ref.  3,  4,  5,  7  in  section  1.3.3). 

Very  few  attention  has  been  given  to  the  heterogeneous  deformation  in  modelling  the 
overall  behavior  of  the  granular  material.  To  account  for  the  heterogeneity,  the  other  three 
methods  listed  above  for  micro-macro  relationship  have  been  developed.  In  the  formulation  of 
static  hypothesis,  the  derived  stress-strain  relationship  is  corresponding  to  a  lower  bound 
solution,  in  contrast  to  the  upper  bound  solution  from  the  kinematic  hypothesis.  Therefore,  the 
development  of  a  static  hypothesis  is  useful  in  providing  a  set  of  bounds  for  estimating  the 
exact  solution.  The  bounds  are  particularly  useful  when  the  microstructure  information  is  • 
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limited  so  that  an  exact  solution  cannot  be  determined. 

To  further  consider  the  heterogeneity,  material  is  divided  into  cells  (i.e.  micro¬ 
elements).  Configurations  of  cells  are  used  to  characterize  the  microstructure.  Self-consistent 
method  was  applied  to  relate  the  micro-element  strain  and  the  macro-scale  strain.  Stress- 
strain  relationship  was  derived  based  on  the  self  consistent  method,  for  a  two  dimensional 
granular  packing  considering  plastic  sliding  (Ref.  3,  8  in  section  1.3.3).  The  method  can  not 
be  readily  extended  to  a  three-dimensional  problem  because  it  requires  a  Green’s  function  for 
general  anisotropic  media  which  is  not  available  for  three-dimensional  condition  in  the  current 
literature.  In  this  project,  an  effort  was  also  devoted  to  the  development  of  general  Green's 
function  (Ref.  11  in  section  1.3.3).  The  developed  three-dimensional  Green's  function  is  a 
useful  fundamental  solution  that  can  be  used  not  only  in  the  self-consistent  method  but  also  in 
other  fields  of  mechanics. 

Statistical  method  was  also  used  to  relate  the  micro-element  strain  and  the  macro¬ 
scale  strain.  The  canonical  ensemble  average  technique  used  in  statistical  mechanics  was 
adopt.  The  problem  was  addressed  and  the  theory  was  derived  considering  particle  sliding 
(Ref.  4,  5,  7  in  section  1.3.3).  Since  this  problem  falls  in  the  category  of  homogenization 
theory,  an  effective  medium  theory  in  this  connection  was  also  examined  (Ref.  21,  30  in 
section  1.3.3). 

The  results  show  that  heterogeneous  strain  field  is  a  very  important  factor.  Micro¬ 
macro  relationship  without  consideration  of  heterogeneous  strain  field  (e.g.  in  the  kinematic 
hypothesis  model)  fails  to  model  yielding  of  a  packing.  The  three  aforementioned  methods  are 
capable  of  modelling  heterogeneous  strain  field  and  yielding  of  a  packing.  However,  the  static 
hypothesis  underestimates  while  the  other  two  methods  overestimate  failure  strength. 

Additional  work  still  needs  to  be  focused  in  this  area.  Other  factors  associated  with  this 
problem  also  need  to  be  considered  such  as  the  method  of  residual  force  redistribution  due  to 
particle  sliding. 

1.2.3.  Glasses  of  Micromechanics  Constitutive  Theory 

A  desirable  approach  is  to  view  the  system  of  discrete  particles  as  a  continuum. 
Regarding  to  this  approach,  the  fundamental  question  is:  what  should  be  the  form  of  the 
equivalent  continuum?  Considering  the  distinctive  two  modes  of  particle  movements, 
displacement  and  rotation,  the  equivalent  continuum  model  for  the  discrete  system  can  be 
represented  by  different  types  of  continue,  for  example,  micro-polar  continuum,  high  gradient 
continuum,  non-local  continuum,  etc.  Differed  from  the  classic  continua  in  solid  mechanics, 
these  continua  display  effects  of  internal  characteristic  length.  When  the  particle  size  is 
relatively  large  compared  to  the  sample  size,  the  notion  of  classic  continua  is  no  longer 
adequate  to  represent  the  granular  system. 
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It  is  clear  that  a  micro-polar  continuum  captures  the  mode  of  particle  rotation,  thus 
models  the  consequence  of  asymmetric  stress  in  granular  material.  It  was  found  that  the 
stress  distribution  and  deformation  in  granular  medium  depend  significantly  on  the  rotation 
stiffness  between  particles  (Ref.  2,  12  in  section  1.3.3). 

The  high  gradient  continuum  depicts  the  effect  of  high  strain  gradients.  The  model 
considers  the  strain  difference  of  neighboring  particles,  thus  can  be  regarded  as  a  non-local 
model  in  a  differential  form.  This  type  of  continuum  is  useful  in  representing  phenomena 
associated  with  internal  characteristic  length  of  the  granular  material  (Ref.  14  in  section  1.3.3). 
For  example,  the 

high  gradient  medium  predicts  the  phenomenon  of  decay  and  dispersion  of  short  length 
waves  (Ref.  15  in  section  1.3.3). 

So  far,  studies  of  the  generalized  continuum  models  show  interesting  singular  behavior 
in  elastic  range.  Plastic  behavior  for  these  continue  has  not  yet  been  studied.  However,  it  is 
expected  that  these  generalized  continue  are  useful  in  modelling  localized  deformation  in 
plastic  region  due  to  microstructure  effect.  Therefore  the  generalized  continue  are  potentially 
useful  models  for  the  plastic  and  after  failure  behavior  which  can  not  be  modelled  by  the 
classic  continuum. 

1.2.4  Contact  Law  of  Particles  with  a  Binder 

Most  studies  available  in  the  literature  are  for  granular  materials  with  frictional  contacts. 
Non-linear  inter-particle  contact  behavior  of  Hertz-Mindlin  type  for  frictional  contact  was  used 
in  the  analysis  for  random  packing  of  multi-sized  granules.  For  particles  with  cemented 
contact  such  as  cemented  sand  have  been  previously  modeled  using  a  adhesive-frictional 
contact  model  (Ref.  27  in  section  1 .3.3).  However,  a  more  realistic  model  for  cemented 
contact  is  to  view  the  particles  connected  with  binders  of  finite  thickness  and  finite  contact 
area.  Along  this  line,  a  theoretical  model  was  developed  for  particles  with  elastic  bonds  and 
visco-elastic  bonds. 

1.3  Summary  of  The  Project 

1.3.1  Award  Information 
Award  Number; 

Amount: 

Period: 

Title: 


AFOSR  F49620-92J-0091 
$  179,964 
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Based  on  Micromechanics 
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CHAPTER2 

Mcrostructuke  Characterizatton 


Microstructure  of  granular  material  has  significant  effects  on  the  macro-scale 
constitutive  behavior.  Method  of  microstructure  characterization  for  granular  material  thus 
plays  an  important  role  for  prediction  of  stress-strain  relationship.  Mathematical  representation 
depends  greatly  on  how  a  basic  element  of  granular  microstructure  is  selected.  Treating  the 
particle  group  as  basic  element,  the  packing  structure  of  material  can  be  represented  by  a  set 
of  particle  group  configurations.  In  this  case,  the  basic  element  consists  of  more  than  two 
particles,  thus  provides  more  information  of  microstructure.  Particle  group  represented  by  a 
Voronoi  cell  has  been  previously  used  in  characterizing  microstructure  and  in  developing 
micromechanical  models  of  granular  materials  (Chang  et.  al.  1992,  Chang  1992,  1993). 

In  this  chapter,  Delaunay  tessellation  is  discussed  for  microstructure  characterization 
of  granular  material.  The  distinctive  geometric  feature  for  Delaunay  cells  is  the  tetrahedral 
shape,  all  cells  composed  of  four  faces.  As  oppose  to  the  polyhedral  Voronoi  cells, 
composed  of  different  numbers  of  faces,  the  Delaunay  cells  potentially  are  simpler  to  be 
modelled  and  easier  to  be  represented  statistically.  Furthermore,  the  movement  of  four 
particles  can  be  exactly  represented  by  linear  displacement  and  rotation  fields.  This  fact 
makes  the  linear  fields  good  representations  for  each  tetrahedral  Delaunay  cell  which  is 
associated  with  only  four  particles.  On  the  other  hand,  each  Voronoi  cell  is  usually  associated 
with  more  than  four  particles.  Thus  artificial  constraints  are  involved  when  linear  fields  are 
used  for  representing  the  movement  of  particles. 

We  seek  for  deriving  a  stress-strain  relationship  of  a  representative  volume  of  granular 
material  directly  from  the  behavior  of  two-particle  interaction.  In  order  to  account  the 
microstructure  represented  by  a  group  of  particles,  we  view  the  packing  structure  of  granular 
materials  at  three  different  levels,  namely,  1)  inter-particle  contact,  2)  micro-element,  and  3) 
representative  volume.  Each  pair  of  particles  is  regarded  as  the  basic  unit  of  granular 
material.  A  micro-element,  which  is  associated  with  several  particles,  is  an  elementary  unit  at 
micro-scale  level.  At  micro-element  level,  the  discrete  particle  system  is  transformed  into  an 
equivalent  continuum  system  and  its  behavior  is  described  by  the  continuum  concepts  of 
stress  and  strain.  A  representative  volume  is  defined  as  an  assembly  which  contains  a  large 
number  of  particles  to  be  representative  of  the  granular  material.  A  two  dimensional 
schematic  representation  of  the  three  levels  of  granular  material  is  shown  in  Fig.  2.1. 

In  this  chapter,  we  derive  the  stress-strain  relationship  for  granular  material  based  on 
micromechanics  approach  using  Delaunay  tessellation  method.  The  derived  stress-strain 
relationship  is  illustrated  by  comparing  the  predicted  micro  and  macro  behavior  with  oomputer 

simulation  results.  The  predicted  results  are  also  compared  for  both  Voronoi  tessellation  and 

* 

Delaunay  tessellation  methods. 
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2.1  Discrete  Behavior  at  Contact  Level 

2.1.1  Kinematics  of  Two  Particles 


Particles  in  this  model  are  assumed  to  be  stiff  such  that  the  contact  deformation  are 
small  compared  to  the  size  of  particles.  For  practical  purpose,  the  movement  of  particles  can 
be  described  by  the  kinematics  of  rigid  particles.  There  are  two  modes  of  movement  for  a 

particle:  translation,  and  rotation,  Aoj".  The  superscript  'rf  refers  to  the  n-th  particle. 

Based  on  the  kinem.atics  of  two  rigid  particles  of  convex  shape  as  shown  in  Fig.  2.2,  the 

relative  displacement  A 5™  and  the  relative  rotation  A6^  of  particle  ’m'  to  particle  'li  at  the 

contact  point  'c'  are  given  by 

Asr  =  A«,"-A«,'^-ej,.j(rrAa.;-rrAcO;”)  (2-1) 

Aer  =  A<or-Au?  (2-2) 

nc 

where  is  the  position  vector  measured  from  the  center  of  the  n-th  particle  to  the  contact 


point  'c'  and  the  quantity  is  the  permutation  symbols  used  in  tensor  representation  for 

cross  product  of  vectors.  Combining  tensor  and  matrix  notation,  Eqs.  2.1  and  2.2  can  be 
written  in  a  compact  form  as  follows: 

{Afl,""!  =  [fiTl  { Al//}  -  {AUj")  (2.2 

where  [Dn  is  the  generalized  inter-particle  displacement  vector  and  [u;]  is  the 
generalized  particle-displacement  vector.  (A-n.  iuf]  ,  and  IRy"]  are  defined  as  follows: 

{/),")  = .  ■:  {ir/i  =  J  ;[«,•!=  '  .  <2.4 

er  I®  J 


where  6j^  is  the  Kronecker  delta. 


2.1.2  Inter-particle  Law 

Due  to  movement  of  particles,  forces  and  moments  transmit  through  at  the  inter¬ 
particle  contacts.  The  inter-particle  displacement,  ASf”,  and  the  contact  force,  Af^,  at  the 
inter-particle  contact  can  be  expressed  in  a  general  relationship  of  incremental  form  as 
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(2.5) 


A  A  ft  Wf 

A/j  =  %  A5y 


The  contact  stiffness  tensor  is  given  by 

tr  "Cnrnr^c&r^r^'ro 

where  n,  s  and  f  are  the  basic  unit  vectors  of  the  local  coordinate  system  construct  at  each 
contact.  The  vector  n  is  the  outward  normal  to  the  contact  plane.  The  other  two  orthogonal 

vectors  s  and  t  are  on  the  contact  plane.  The  scaler  quantities  is  the  inter-particle 

normal  stiffness  and  is  the  inter-particle  shear  stiffness. 

The  inter-particle  rolling  and  torsional  stiffness  transmit  the  inter-particle  moment 
(couple).  The  inter-particle  torsional  stiffness  is  denoted  by  and  the  inter-particle  rolling 

stiffness  is  denoted  by  g”'” .  The  inter-particle  relative  rotation,  A0"”*,  and  the  inter-particle 

contact  moment,  A/tt,-"**,  can  be  expressed  in  a  general  relationship  of  incremental  form  as 

(2.7) 


Amr  =  giP  Ae^ 

The  inter-particle  rotation  stiffness  tensor  is  given  by 


wn  yimnmnm  wn/  /un  nm 
8ij  =  8n  «y  is-i  Sj  Hi  tj  ^ 

Eqs.  2.5  and  2.7  can  be  combined  into  a  compact  form  as  follows: 


(2.8) 

(2.9) 


{A^i*™}  =  [Kf] 

where  [Dn  is  the  generalized  inter-particle  displacement  vector,  {Fn  is  the  generalized 
inter-particle  force  vector,  and  IKP]  is  the  generalized  inter-particle  stiffness  tensor,  defining 


in  the  follows: 


{D^}  = 


l/r] 

% 

0  ' 

■;  IKf]  = 

V 

arm 

j 

r‘ 

0 

_  ipn 

(2.10) 


2.2  METHODS  OF  MiCROSTRUCTURE  CHARACTERIZATION 


To  study  a  granular  material  with  randomly  packed  particles,  a  suitable  model  for 
characterizing  microstructure  is  essential.  Two  simple  and  widely  used  methods  in  the  theory 
of  stochastic  geometry  are  examined,  namely,  the  Voronoi  tessellation  method  and  the 
Delaunay  tessellation  method. 
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Both  methods  was  originally  used  to  characterize  the  microstructure  of  a  set  of 
randomly  distributed  points  in  space,  such  as  the  two  dimensional  example  schematically 
shown  in  Fig.  2.3a.  The  configuration  of  a  set  of  discrete  points  has  certain  statistical  similarity 
to  that  of  a  set  of  discrete  particles,  such  as  the  two  dimensional  example  schematically 
shown  in  Fig.  2.3b.  Therefore,  the  methods  for  a  set  of  points  are  expected  to  have  ability  for 
capturing,  not  in  detail  but,  main  features  of  the  microstructure  for  a  set  of  randomly 
distributed  particles  in  space. 

Both  methods  transform  the  set  of  points  into  a  cell  structure.  Microstructure  of  the 
system  is  characterized  through  the  arrangements  and  shapes  of  cells.  The  Voronoi 
tessellation  method  divides  the  space  into  polyhedral  cells.  Each  Voronoi  cel!  contains  a 
particle  point.  Each  face  of  the  polyhedron  separates  two  neighboring  particles.  Shape  of  the 
polyhedron  represents  the  arrangement  of  a  particle  point  and  its  neighbors.  The  Voronoi 
tessellation  for  the  set  of  points  in  Eig.  2.3a  is  given  in  Fig.  2.4a.  If  each  polyhedron  cell  is 
occupied  by  a  particle,  it  is  corresponding  to  Fig.  2.3b. 

The  Delaunay  tessellation  is  a  unique  representation  dual  to  the  Voronoi  tessellation. 
The  Delaunay  method  divides  the  space  into  tetrahedral  cells.  Based  on  the  vertexes  of  the 
Voronoi  polyhedra,  the  Delaunay  cell  can  be  formed  by  connecting  the  four  particle  points 
nearest  to  the  vertex.  Each  tetrahedral  cell  has  4  faces  and  6  edges.  The  cell  encloses  a 
void.  Each  edge  connects  a  pair  of  particles.  The  Delaunay  tessellation  for  the  set  of  particle 
points  in  Fig.  2.3a  is  given  in  Fig.  2.4b. 

The  method  of  Voronoi  tessellation  has  been  utilized  by  Chang  et.  al.  (1992)  and 
Chang  (1992,  1993)  to  characterize  microstructure  of  granular  material.  The  Delaunay 
tessellation  is  expected  to  have  advantages  over  the  Voronoi  tessellation  in  that  all  Delaunay 
cells  are  tetrahedra,  of  which  the  statistical  information  are  easier  to  be  obtained.  In  this 
chapter,  we  propose  the  Delaunay  tessellation  for  characterizing  microstructure. 


2.3  Delaunay  Tessellation  for  Micro-element 


2.3.1  Strain  Field  in  a  Delaunay  Micro-element 


A  Delaunay  cell  is  associated  with  four  particles.  The  tetrahedral  micro-element  is 
constructed  by  the  centroids  of  the  four  particles.  Define  displacement  and  rotation  fields  in 
the  micro-element  as  linear  functions; 


LufX)  =  (2.11) 

A<Oj(J£)  =  Ato^+Aoi^X^.  (2.12) 
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where  Xj  is  the  position  veotor  measured  from  the  centroid  of  the  micro-element,  Ui  and  oi- 

are  respectively  the  displacement  and  the  rotation  at  the  centroid  of  the  micro-element,  is 

the  displacement  gradient,  <0*^-  is  the  rotation  gradient,  and  the  superscripts  'e'  refers  to  the  e- 
th  micro-element. 

Using  continuum  variables  of  the  displacement  and  rotation  fields,  the  inter-particle 
displacement  and  the  inter-particle  rotation  in  Eqs.  2.1  and  2.2  can  be  expressed  as: 

A6r  =  P.13) 

Aer  = 

where  l””*  is  the  branch  vector  measured  from  the  center  of  the  n-th  particle  to  the  center  of 

the  m-th  particle  and  is  the  position  vector  connecting  the  center  of  n-th  particle  to  the 

contact  point  'o’.  The  branch  vector  and  position  vectors  are  related  by  or 

=  rj^-rJ^ .  The  asymmetric  deformation  strain  and  the  polar  strain  are  defined, 
similar  to  that  given  by  Chang  and  Liao  (1990)  and  Chang  and  Ma  (1990),  as  follows: 

=  Ax,.'.  +  AuJ  (2' '5) 

Ay^  =  A<o',- 

The  symmetrical  part  of  the  strain  Ae^-  is  equal  to  the  symmetrical  part  of 
displacement  gradient,  representing  the  usual  strain  of  the  micro-element,  i.e.. 

The  skew  symmetric  part  of  the  strain  Ae|  represents  the  net  spin  of  particles  (i.e.,  the 

difference  between  rigid  body  rotation  and  the  particle  rotation  at  the  centroid  of  the  micro¬ 
element).  Thus 

^4]  =  \  =  I  ^ 

It  is  noted  that,  we  aim  to  describe  the  movements  of  discrete  particles  with  continuum 
variables  {  4.  y;-}  instead  of  discrete  variables  {«A  to”} .  Substituting  the  linear  fields  of 


12 


Eqs.  2.15  and  2.16  into  Eqs.  2.11  and  2.12,  particle  movements  and  micro-element  strains  are 
related  by 


A«,(X“)  =  A«,%Aef.X,"-ej,.^Qtx; 

(2.19) 

Au^X”)  =  Acof  +  AT^X/ 

(2.20) 

where  W,”,  co"  are  respectively  the  displacement  and  rotation  of  particle  n.  X”  is  the  position 

vector  from  centroid  of  the  micro-element  to  the  centroid  of  particle  n.  Since  the  superscript  n 
denotes  the  particles  from  1  to  4  for  a  Delaunay  cell,  the  set  of  24  discrete  variables 

can  be  fully  replaced  by  the  set  of  24  continuum  variables  «f} 


Neglecting  the  rigid  body  translation  and  the  rigid  body  rotation,  we  let 

Aut  =  0  (2-21) 

AUij-LUj^  =  0  (2.22) 

Substituting  Eq.  2.22  into  Eq.  2.18,  it  becomes: 

Aaf  -  (2.23) 

Using  Eqs.  2.21  and  2.23,  Eqs.  2.19  and  2.20  can  now  be  written  as  follows; 

Ab/J:")  =  i(A4<-Ae;iX;  (2.24) 

Au,(X')  -  (2.25) 


In  Eqs.  2.24  and  2.25,  the  total  number  of  continuum  variables  is  18  as  represented  by 
the  set  {4  41  With  those  6  known  conditions  of  rigid  body  movement  (Eqs.  2.21  and 

2.23),  the  24  discrete  variables  of  can  be  uniquely  defined  from  the  set  of 

18  continuum  variables  {e-,  y’-} 

It  is  noted  that  the  linear  field  assumption  used  in  Voronoi  representation  (Chang, 

1992)  does  not  lead  to  a  unique  mapping  between  discrete  variables  and  continuum 
variables.  As  oppose  to  the  Delaunay  cell,  a  Voronoi  cell  usually  consists  of  more  than  4 

particles,  such  that  the  number  of  discrete  variables  are  more  than  24. 

Therefore,  the  linear  field  assumption  in  Voronoi  representation  involves  artificial  kinematic 
constraints.  In  this  regards,  linear  field  assumption  is  a  better  representation  for  Delaunay 
tessellation  method. 
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2.3.2  Micro-element  Strain  versus  Inter-particle  Displacement 


We  express  the  relative  movement  of  a  particle  pair  (Eqs.  2,13  and  2.14)  in  terms  of 
the  generalized  strain  tensor  {Efk)  as  follows: 

{Acri  =  {AE^l  P.26) 

where  {Dr)  ,  as  given  in  Eq.  2.10,  is  the  generalized  inter-particle  displacement  vector 


between  particles 'm'  and  'n'  and  \-L^]  is  the  fabric  tensor  operator.  and  are 

given  by 


V;  -  /vtn  me 

;  [Ejl)  =  < 

0  b.jr 

e 

M 

(2.27) 


Each  Delaunay  cell  of  tetrahedral  shape  consists  of  4  vertices  and  6  branch  vectors 
(edges).  Since  there  are  6  variables  of  associated  with  each  branch  vector,  thus  the 


total  number  of  variables  {A"”}  is  36  for  a  cell,  which  can  be  determined  from  18  continuum 


variables  of  thfe  set  {e^.,  •  It  is  noted  that  the  36  variables  are  not  independent  due  to 

geometric  constraint  of  compatibility.  On  the  face  made  of  3  vertices  (say,  m,  n,  and  q),  the  6 
variables  of  [Dn  associated  with  the  branch  vector  can  be  expressed  by  the  variables 

and  [Dn  associated  with  the  other  two  branch  vectors  and  given  by 

Dr  =  nr+Dr 

For  each  Delaunay  cell,  18  compatibility  equations  can  be  written.  The  other  12  equations 
are: 


A**  =  a'^+A^ 


(2.29) 


The  inter-particle  displacements  obtained  from  the  generalized  strain  tensor  {Eft)  using  Eq. 
2.26  automatically  satisfy  the  18  compatibility  equations. 


2.3.3  Inter-particle  force,  versus  micro-element  stress 

The  stress  of  a  micro-element  cari  be  defined  using  the  principle  of  energy 
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equivalence.  Due  to  an  increment  of  strain,  the  potential  energy  from  the  stress  and  strain  is 
equal  to  the  work  done  at  the  inter-particle  contacts.  Note  that  each  particle  pair  is  the  edge 


shared  by  several  micro-elements.  Let  tx''  be  the  number  of  micro-elements  that  share  the 
same  edge.  It  is  noted  that,  for  simplicity,  we  use  the  superscript  'c'  to  instead  of  the 
superscripts  'nm’.  Assuming  the  work  done  at  this  contact  is  equally  shared  by  all  micro¬ 
elements  associated  with  this  contact.  The  energy  equivalence  is 


OjAeJ  turA-i-r  "  iSf+m/Aef) 


(2.30) 


ce¬ 

lt  can  be  written  in  a  compact  form  as  follows; 


{S;f{AE^]  =  —  E  —  (2-31) 

V*  c 


where  is  the  generalized  stress  tensor  given  by 


(2.32) 


in  which  {o^}  is  the  Cauchy  stress  tensor  and  is  the  polar  stress  tensor  for  a  micro¬ 


element  It  is  noted  that  is  2  for  two  dimensional  case.  In  a  three  dimensional  case,  is 
an  additional  measure  of  microstructure  describing  the  neighbors  of  a  micro-element. 

By  substituting  Eq.  2.26  into  Eq.  2.31, 


{AB'l 

I  yc 


cur 


=  0 


(2.33) 


Since  the  incremental  strain  is  arbitrary  chosen,  in  order  to  satisfy  the  energy 
conservation,  we  have 


c  (xr 

(2.34) 

c  a 

(2.35) 

More  specifically. 
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(2.36) 


F"  T 


(2.37) 


c  a" 


where  -  l^  -  X^-Xl*.  Also,  fj  -  and  mf  =  nij”*  are  the  inter-particle  force 
and  moment  acting  on  the  r>th  particle  respectively. 


2.3.4  Constitutive  Law  for  Micro-element 


The  stress-strain  relationship  for  a  micro-element  is  given  by 

(ASj)  =  [CjJ  {hE’u)  (2-38) 

The  stiffness  tensor  can  be  derived  from  the  following  three  essential  relationships: 

1)  inter-particle  force  to  micro-element  stress  relationship  (Eq.  2.35) 

2)  inter-particle  law  (Eq.  2.9) 

3)  micro-element  strain  to  inter-particle  displacement  relationship  (Eq.  2.26) 

These  relationships  are  schematically  shown  in  Fig.  2.5.  Using  these  three  equations,  the 

stiffness  tenson'  is  derived,  given  by 


c 

Eq.  2.38  can  be  expressed  as  follows: 


where 


lAoJ. 


^ijki 


Ae 


tf 


c 

_  1  1  _  fCjL  /V^  ^ 

c 

D’m  =  ^E^w.vxrc-*”'-."') 


C  a'^ 

V  c  Ct 

A  brief  summary  of  the  discrete  variables,  the  continuum  variables,  and  their 
relationships  for  a  Delaunay  micro-element  is  given  in  Table  2.1. 


(2.39) 


(2.40) 


(2.41) 
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2.4  Overall  Behavior  for  Representative  Volume 

In  what  follows,  we  derive  the  constitutive  law  for  a  representative  volume  of  the 
granular  material.  The  macro-level  stress-strain  behavior  is  obtained  from  the  averaged 

behavior  of  micro-elements.  The  overall  stress  and  strain,  denoted  by  AS^j  and  ,  are 

regard  as  volume  averages  of  the  local  stresses  and  local  strains  at  the  micro-element  level, 
such  that: 

(2.42) 

^  e 

-  4E  (2.43) 

^  e 

where  F*  is  the  volume  associated  with  the  e-th  micro-element.  Summation  over  the  volume 

F®  of  all  micro-elements  is  equal  to  the  representative  volume  F. 

To  derive  the  correct  overall  constitutive  relationship,  it  is  essential  to  account  for  the 

fluctuations  of  local  strain  A£"  and  local  stress  as;  of  micro-elements  in  a  heterogeneous 

granular  system.  The  fluctuations  are  complicated  in  nature  and  difficult  to  model.  Strain 
fluctuations  for,  granular  material  without  consideration  of  polar  strain  has  been  studied  using 
Self-consistent  method  (Chang  et.al.  1992)  and  using  method  of  Distributive  tensor  (Chang 
1993).  It  is  noted  that  we  will  evaluate  the  method  of  Distributive  tensor  in  detail  for  predicting 
the  mechanical  behavior  of  heterogeneous  granular  material  in  Chapter  3.  In  trade  off 
complicate  analyses,  it  is  desirable  to  conduct  simple  analyses  for  obtaining  approximate 
solutions.  In  some  situations,  the  results  of  simple  analyses  involve  straight  forward  and 
tractable  assumptions  which  may  be  more  useful  in  providing  insights  on  the  nature  of 
fluctuation  of  local  stress  and  strain.  Although  the  results  of  simple  analyses  are  approximate 
to  the  exact  solutions,  they  provide  bounds  and  estimates  which  are  useful  for  assessing  the 
range  of  mechanical  behavior  for  the  complex  granular  systems.  To  this  end,  we  conduct 
analyses  with  the  following  four  simple  assumptions:  (For  clarity  in  the  subsequent 
presentation,  we  drop  the  subscript  for  expressions  of  stress,  strain  and  stiffness  tensors.) 

(1)  uniform  strain  (i.e.,  €,Y  aie  constants  for  all  micro-elements) 

(2)  uniform  stress  (i.e.,  o,  p.  are  constants  for  all  micro-elements) 

(3)  uniform  strain  and  uniform  polar  stress  (i.e.,  e,  (1  are  constants  for  all  micro¬ 
elements) 

(4)  uniform  stress  and  uniform  polar  strain  (i.e.,  o,  Y  are  constants  for  all  micro¬ 
elements) 
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Clearly,  the  uniform  strain  assumption  of  case  (1)  leads  to  an  upper  bound  solution, 
the  uniform  stress  assumption  of  case  (2)  leads  to  a  lower  bound  solution,  the  partial  uniform 
strain  or  stress  assumptions  (i.e.,  cases  (3)  and  (4))  lead  to  intermediate  solutions.  The  overall 
stiffness  tensor  can  then  be  derived  based  on  the  four  assumptions,  as  follows: 

Case  (1)  -  uniform  strain  assumption  (€,  y  are  constants);  The  effective  macro-scale 
stiffness  tensor  is  the  volume  average  of  stiffness  tensor  for  all  micro-elements,  given  by 

=c  = (2.44) 

^  e 

where  the  overline  is  defined  as  volume  average  of  local  quantities.  This  definition  of  overline 
applies  throughout  this  section. 

Case  (2)  -  uniform  stress  assumption  (o,  jl  are  constants):  Inverse  of  the  effective 
macro-scale  stiffness  tensor  is  the  volume  average  of  the  inverse  of  stiffness  tensor  for  all 
micro-elements,  given  by: 

(2.45) 

where  the  overline  of  C~^  is  defined  as  the  volume  average  of  the  inversed  local  stiffness. 


Case  (3)  -  uniform  strain  and  uniform  polar  stress  assumption  (e,  p.  are  constants): 
The  effective  stiffness  tensor  is  given  by: 


= 


(A-^y^  (A-^yu-^B 

DA  ~^(A^y^A^  +E-DA 


(2.46) 


where  the  overline  of  is  defined  as  the  volume  average  of  the  inversed  local  stiffness  A^. 
The  expressions  of  A®,  5*^,  and  £*  can  be  found  in  Eq.  2.41. 

Case  (4)  -  uniform  stress  and  uniform  polar  strain  assumption  ( O ,  Y  are  constants): 
The  effective  stiffness  tensor  is  given  by: 


= 


(E^y^E-^D 


BE-^(E-^y^ 

(E^r 


(2.47) 


2.5  Model  Performance 


To  evaluate  the  applicability  of  the  derived  micromechanical  model  with  Delaunay 
tessellation  characterization,  we  study  elastic  stress-strain  behavior  of  a  random  packing  of 
planar  disks.  The  predicted  results  are  compared  with  the  computer  simulation  results. 

A  periodic  space  of  randomly  packed  particles  used  in  this  study  is  shown  in  Fig.  2'.6. 
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The  periodic  packing  can  be  repetitively  stacked  to  fill  the  entire  space,  thus  is  a 
representative  volume  of  the  granular  material.  The  packing  is  composed  of  1200  particles  of 
radius  0.175  mm  and  800  particles  of  radius  0.100  mm.  The  total  number  of  contacts  in  the 
packing  is  5330.  The  coordination  number  (the  average  number  of  contacts  per  particle)  is 
5.33,  the  area  of  the  representative  volume  is  174  mm^  and  the  void  ratio  is  0.234. 

Mechanical  behavior  of  the  packing  is  obtained  from  a  computer  simulation  method 
using  the  discrete  element  method  by  Chang  and  Misra  (1989a).  The  packing  is  divided  into 
2000  Voronoi  cells  shown  in  Fig.  2.7  and  divided  into  4070  Delaunay  cells  shown  in  Fig.  2.8. 
The  predicted  results  obtained  from  the  present  model  using  Delaunay  tessellation  method  are 
first  compared  with  the  computer  simulation  results  and  then  compared  with  the  results 
obtained  from  Voronoi  tessellation  method. 

2.5, 1  Comparison  witli  Computer  Simulation 

To  illustrate  the  heterogeneous  nature  of  strain  and  stress  in  the  granular  system,  the 
variation  of  these  quantities  are  investigated  at  the  local  level.  The  loading  condition  is  taken 
to  be  pure  shear.  The  overall  stress  and  strain  for  the  representative  volume  are  specified  to 
be:  s„  =  0.22%,  =  -0.22%,  -  0,  Gi,2,  =  0,  p,3  =  0,  and  =  0.  The  overall  stress  and 

strain  for  the  representative  volume  to  be  computed  are:  a,,,  cr^^,  Ec^;,  E^,^;,  y,^,  and  723.  Also 
computed  in  the  model  are  the  local  stress  and  local  strain  for  each  micro-element.  The  inter¬ 
particle  stiffness  are  assumed  to  be  same  for  all  inter-particle  contacts,  k„  =  50  kN/m,  =  40 
kN/m,  and  =  0.  The  material  of  concerned  is  composed  of  round  particles  with  no  rotational 
resistance. 

We  list  the  overall  stresses  and  strains  for  the  representative  volume  in  Table  2.2  and 
the  standard  deviation  of  local  stresses  and  local  strains  for  micro-elements  in  Table  2.3, 
predicted  under  the  four  different  assumptions  of  heterogeneous  states.  The  results  from 
computer  simulation  are  also  listed  for  comparison.  According  to  Tables  2.2  and  2.3,  the 
results  based  on  the  assumptions  of  uniform  polar  strain  (cases  1  and  3)  are  in  good 
agreement  with  those  from  the  computer  simulation.  However,  the  results  based  on  the 
assumptions  of  uniform  polar  stress  (cases  2  and  4)  deviate  greatly  with  those  from  the 
computer  simulation.  Therefore,  the  assumption  of  uniform  polar  strain  seems  to  be  better  for 
granular  material. 

The  predicted  effective  moduli  are  compared  with  those  obtained  from  computer 
simulation  to  evaluate  the  micromechanical  method.  Figures  2.9  and  2.10  show  the  shear  and 
bulk  moduli  with  different  ratio  of  inter-particle  stiffness  kjk„,  under  four  different  assumptions 
of  the  heterogeneous  state  of  the  granular  packing.  In  this  prediction,  for  all  inter-particle 
contacts,  the  normal  inter-particle  stiffness  /f„  =  50  kN/m,  while  tangential  inter-particle  stiffness 
k^  is  chosen  as  0,  10,  20,  30,  40,  and  50  kN/m.  In  both  Figures,  case  1  is  the  upper  bound" 
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solution  and  case  4  is  the  lower  bound  solution.  Case  2  with  uniform  polar  stress  does  not 
give  good  agreement,  in  which,  the  value  of  has  little  effect  on  moduli.  The  results  based 
on  the  assumptions  of  uniform  polar  strain  (cases  1  and  3)  provide  a  good  range  for  upper 
and  lower  estimates  respectively. 

2.5.2  Comparison  with  Voronoi  Tessellation  Method 

The  predicted  effective  moduli  based  on  Delaunay  tessellation  method  are  compared 
with  those  based  on  Voronoi  tessellation  method.  Figures  2.11  and  2.12  show  the  shear  and 
bulk  moduli  with  different  ratio  of  inter-particle  stiffness  kjk„.  In  both  figures,  the  predicted 
moduli  based  on  Delaunay  tessellation  method  are  identical  with  those  predicted  based  on 
Voronoi  tessellation  method  for  case  1  while  the  predicted  moduli  based  on  Delaunay 
tessellation  method  are  lower  than  those  predicted  based  on  Voronoi  tessellation  method  for 
case  3.  This  is  due  to  the  artificial  constraint  involved  in  the  kinematic  assumption  in  Voronoi 
tessellation  method  as  described  previously.  The  predicted  moduli  from  case  1  and  case  3 
for  Delaunay  tessellation  method  present  a  range  to  cover  the  computer  simulation  results. 

2.6  Summary 

In  the  rnicromechanics  approach,  structure  of  the  granular  material  can  be 
represented  by  a  set  of  micro-elements  characterized  by  Delaunay  or  Voronoi  tessellation.  We 
transform  the  discrete  variables  of  particle  translation  and  particle  rotation  into  continuum 
variables  of  strain  and  polar  strain  of  micro-elements.  Correspondingly,  we  transform  the 
discrete  variables  of  inter-particle  force  and  moment  into  continuum  variables  of  stress  and 
polar  stress. 

In  a  tetrahedral  Delaunay  cell,  the  assumption  of  linear  displacement  and  rotation 
fields  leads  to  a  unique  mapping  between  discrete  variables  and  continuum  variables.  This 
uniqueness  is  not  guaranteed  for  a  polyhedral  Voronoi  cell  comprised  of  more  than  four  faces. 
Thus  linear  field  assumption  used  in  Voronoi  tessellation  induces  artificial  constraints. 

The  macro-scale  constitutive  relationship  for  a  representative  volume  can  be  derived 
using  four  simple  assumptions  on  the  fluctuation  of  micro-element  stress  and  micro-element 
strain.  Based  on  the  example  packing,  the  assumption  of  uniform  polar  strain  of  the  present 
model  gives  results  comparable  with  those  obtained  from  computer  simulation  method. 
However,  uniform  polar  stress  seems  not  to  be  a  good  assumption.  With  the  assumptions  of 
uniform  strain  and  uniform  polar  strain,  the  present  model  provides  an  upper  bound  solution. 
With  the  assumptions  of  uniform  stress  and  uniform  polar  strain,  the  present  model  provides  a 
lower  estimate  solution. 
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Discrete  variables; 

Generalized  particle  displacement 

{U-] 

Translation  {«,"} 

Rotation  {co"} 

Generalized  inter-particle  displacement 

(A"} 

Inter-particle  displacement  {6-} 

Inter-particle  rotation  {0;} 

Generalized  inter-particle  force 

Inter-particle  force  {fi) 

Inter-particle  moment  {rtii] 

Inter-particle  constitutive  equation  =  [K^^{Dj} 

Continuum  variables: 

Generalized  strain  {F,y} 

Stretch  strain  {cy} 

Polar  strain  {y^} 

Generalized  stress  {S,y} 

Cauchy  stress  { Oy} 

Polar  stress  {py} 

Constitutive  equation  {S,y}  =  [C,yy]{Fy} 

Note: 

_0  Sy 

c 

V/  e^iX-rr-X^r-) 

0 

Table  2.1  Summary  of  the  discrete  variables,  the  continuum  variables,  and  their 
relationships  for  a  Delaunay  micro-element  - 
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Table  2.3  Predicted  standard  deviations  of  local  stresses  and  local  strains  for 
Delaunay  micro-elements 
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inter-particle  contact 

REPRESENTATIVE  VOLUME 


Figure  2.1  Schematic  figure  for  three  levels  of  granular  materials 
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(b) 


Figure  2.4  Microstructure  characterization;  (a)  Voronoi  tessellation;  (b)  Delaunay 
tessellation 
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Micro-element 
constitutive  law 


Micro-element 
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Micro-element 
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Inter-particle 
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Figure  2.5  Relationships  between  kinematic  and  static  variables  for  a  Delaunay 
micro-element 
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Figure  2.6  Packing  structure  of  the  disk  assembly  used  in  example 
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Figure  2.7  Voronoi  tessellation  of  the  disk  assembly 
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Figure  2.8  Delaunay  tessellation  of  the  disk  assembly 
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Figure  2.9  Predicted  effective  shear  moduli  using  the  Delaunay  tessellation  method 
under  four  different  assumptions 
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Figure  2.10  Predicted  effective  bulk  moduli  using  the  Delaunay  tessellation  method 
under  four  different  assumptions 
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Figure  2.11  Comparison  of  effective  shear  moduli  predicted  by  methods  based  on 
Voronoi  tessellation  and  Delaunay  tessellation 
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Bulk  modulus,  B  (kPa) 


Figure  2.12  Comparison  of  effective  bulk  moduli  predicted  by  methods  based  on 
Voronoi  tessellation  and  Delaunay  tessellation 
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CEIAFIER3 

Classes  of  Eoiwaient  Conitnua  for 
Granular  Materials 


It  is  desirable  to  describe  the  deformation  behavior  of  granular  materials  using 
continuum  variables  such  as  strain  and  stress.  However  the  micro  scale  deformation, 
attributed  to  the  relative  movements  of  discrete  particles,  is  not  a  continuous  field.  Therefore, 
for  the  purpose  of  analysis,  it  is  a  fundamental  step  to  transform  the  discrete  particle  system 
to  an  equivalent  continuum  system.  Choice  of  the  transformation  methods  results  to  different 
types  of  equivalent  of  continuum.  This  chapter  describes  different  classes  of  continuum 
which  are  suitable  for  representing  granular  media.  More  complex  continue  capture  more 
features  of  the  deformation  behavior  for  granular  material.  The  continue  of  micropolar  type 
reflect  the  mechanism  of  moment  transmitting  in  granular  media  while  the  high-gradient 
continua  capture  the  effect  of  particle  size  on  wave  dispersion. 


3.1  Introduction 

Granular  material,  such  as  soil,  powder,  ceramic  material,  etc.,  can  be  perceived  as  a 
collection  of  particles.  The  overall  mechanical  properties  for  granular  materials  depend 
significantly  on  the  micro-scale  geometric  arrangement  and  on  the  contact  stiffness  between 
two  interacting  particles.  Therefore  the  effects  of  particle  interaction  should  be  considered  in 
the  construction  of  a  mathematical  model.  However,  a  granular  system  is  so  intricate  in  its 
details  so  as  to  become  essentially  intractable.  Predictions  of  each  particle  movement  are 
costly  to  extract  and  usually  contain  information  well  in  excess  of  one’s  needs.  It  is  therefore 
desirable  to  transform  a  discrete  particle  systems  into  an  equivalent  continuum  system.  The 
deformation  behavior  of  the  granular  material  is  described  using  suitable  continuum  variables 
instead  of  the  movement  of  individual  particles.  Granted,  the  continuum  variables  can  be 
used  to  retrieve  and  keep  track  the  movement  of  each  discrete  particles,  at  least  in  an 
approximate  manner.  This  approach  allows  the  construction  of  a  macro-scale  constitutive 
relations  for  the  equivalent  continuum  of  granular  materials.  The  constitutive  relations  of 
continuum  type  is  capable  of  reflecting  the  discrete  nature  of  particle  interaction  at  micro¬ 
scale  level. 

A  number  of  studies  have  been  attempted  along  this  line  of  approach.  For  example, 
stress-strain  relationship  was  modelled  for  regularly  packed  elastic  spheres  by  Deresiewicz 
(1958),  Duffy  (1959),  Duffy  and  Mindlin  (1957),  and  more  recently,  for  randomly*packed 
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granules,  by  Koenders  (1987,  1994),  Cambou  (1993),  Bathurst  and  Rothenberg  (1988), 

Chang  (1988),  Jenkins  (1988,  1993),  Walton  (1987),  etc.  These  stress-strain  models  treat  the 
granular  material  as  the  traditional  continuum  in  solid  mechanics.  The  stress-strain 
relationship  defined  for  an  infenitesimal  element  in  solid  mechanics  is  now  defined  for  a 
representative  unit  cell  of  the  granular  media.  Since  the  representative  unit  cell  must  contain 
a  large  number  of  particles  to  be  representative  of  the  material,  the  cell  is  of  finite  volume  and 
the  stress  and  strain  actually  fluctuate  within  the  representative  unit  cell.  Thus  the  so  called 
stress-strain  relationship  are  referred  to  the  average  stress  and  average  strain  of  the  unit  cell. 

To  describe  in  more  detail  on  the  deformation  of  a  representative  unit,  it  is  necessary  to 
include  suitable  continuum  descriptors  other  than  stress  and  strain.  For  example,  with  the 
consideration  of  particle  rotation  gradients,  the  equivalent  continuum  resembles  the  continua 
of  micro-polar  or  Cosserat  type  (Chang  and  Liao  1990,  Chang  and  Ma  1991,  1992). 
Corresponding  to  the  rotation  gradients,  moments  are  capable  to  transmit  through  the 
granular  medium,  which  represents  a  mechanism  that  can  not  be  modelled  by  the  traditional 
simple  continuum.  The  effect  of  moment  transmitting  becomes  important  for  particle 
assemblies  with  cemented  inter-particle  contacts  or  with  angular  particles  where  the 
magnitudes  of  moment  transmitted  through  contacts  are  significant  (Chang  and  Ma  1992). 

With  the  consideration  of  higher  orders  of  displacement  gradients,  the  equivalent 
continuum  resembles  the  continua  of  high-gradient  or  non-local  type.  The  use  of  higher 

f 

orders  of  displacement  gradients  reflects  the  effect  of  particle  size  on  the  propagating 
velocities  of  waves  through  granular  media  (Chang  and  Gao  1995b),  which  is  significant  for 
granular  medium  with  large  particle  sizes  (or  shortwave  lengths).  This  dispersion 
phenomenon  of  wave  can  not  be  portrayed  by  the  traditional  simple  continuum.  Elucidated 
from  the  recent  studies  in  high  gradient  models  for  metals  (Coleman  and  Hodgdon  1985, 
Triantafyllidis  and  Aifantis  1986,  Bazant  and  Pijandier-Cabot  1987,  Bardenhagen  and 
Triantafyllidis  1994),  the  model  may  also  be  useful  in  the  study  of  localized  deformation 
associated  with  shear  band  for  granular  media. 

In  this  chapter,  we  first  represent  the  displacement  and  rotation  fields  by  continuous 
polynomial  functions.  Based  on  the  two  continuous  fields,  we  describe  different  classes  of 
equivalent  continua.  To  illustrate  the  relationship  between  micro-macro  properties  such  as 
internal  length  of  granular  material  and  inter-particle  stiffness,  we  discuss  the  constitutive 
constants  of  a  high-gradient  model  for  isotropic  granular  materials.  Example  is  also  given  to 
show  the  links  between  the  high-gradient  model  and  nonlocal  model. 

3.2  Kinematic  Description  of  Granlilar  Material 

A  simple  conceptual  model  for  granular  material  is  to  treat  it  as  a  collection  of  particles. 
When  a  representative  volume  element  of  particles  is  subjected  to  an  increment  of  load. 
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particles  undergo  translations  uf  and  rotations  (Oj  ,  resulting  in  the  relative  deformation 

between  particles.  Based  on  the  kinematics  of  two  particles,  the  inter-particle  displacement 
and  rotation  due  to  relative  movement  of  two  particles  generate  inter-particle  forces  and  inter¬ 
particle  moments  in  granular  materials.  The  magnitudes  of  forces  and  moments  depend  on 
the  contact  stiffnesses. 

For  two  identical  elastic  spherical  particles  in  contact,  the  contact  area  is  of  circular 
shape,  values  of  contact  stiffness  can  be  obtained  from  Hertz-Mindlin  contact  theory  as  a 
function  of  particle  properties  and  contact  force.  While  the  translational  contact  stiffness  are 
responsible  for  transmitting  forces  through  inter-particle  contacts,  the  rotational  contact 
stiffness  are  responsible  for  transmitting  moments  through  inter-particle  contacts  in  a  granular 
media. 

To  develop  a  continuum  mechanics  model  for  the  behavior  of  particle  assembly,  it  is 
desirable  to  view  the  discrete  system  as  an  equivalent  continuum  system.  To  this  end,  we 
define  the  displacement  and  rotation  of  discrete  particles  as  macro-scale  continuum  fields, 

denoted  as;  displacement  u_£  {x)  and  rotation  (x)  .  When  x  is  located  at  the 

centroid  of  n-th  particle,  the  functions  represent  the  displacement  and  rotation  of  the  n-th 
particle. 


(3.1) 


Following  the  approach  by  Chang  and  Liao  (1990)  using  polynomial  expansion  for  the 
displacement  and  rotation  fields,  the  displacement  and  rotation  of  the  m-th  particle  can  be 
estimated  from  the  deformation  gradients  at  the  neighbouring  n-th  particle,  given  by 


,111  HBt  wit  , 

«,  -Ui  *UijLj  ■•■■■ 

Wj  +<i)y£y  +  — WyjLy  Lf.  + 


(3.2) 


where  the  branch  vector  =-r^ ;  the  displacement  and  rotation  derivatives  of  the  n- 

th  particle  are  treated  to  be  the  continuum  kinematic  variables. 

It  is  noted  that  the  rotation  and  displacement  are  not  two  independent  variables. 

During  deformation  of  a  granular  material,  the  rotation  of  a  particle  consists  of  two 
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components:  the  rigid  body  rotation  and  the  particle  spin.  While  the  particle  spin  is 
independent  of  the  displacement  field,  the  rigid  body  rotation  is  solely  induced  by  the 
displacement  field  from  the  skew-symmetric  part  of  displacement  gradients. 

3.3  Different  Glasses  of  Gontinua 

The  discrete  system  can  be  characterized  as  an  equivalent  continuum  system  based 
on  the  two  continuous  fields  of  displacement  and  rotation.  The  discrete  variables  and  the 
continuum  variables  are  linked  through  the  series  expression  of  Eq.  3.2.  When  the  number  of 
terms  in  this  series  is  large  enough,  the  continuum  system  can  capture  as  much  details  as 
the  discrete  system.  A  simplified  kinematic  description  leads  to  a  simpler  continuum  model 
for  the  discrete  system  while  still  capturing  the  essential  features  of  granular  materials. 

Different  degrees  of  approximation  of  Eq.  3.2  lead  to  different  classes  of  ’continua’.  We 
classify  the  models  into  categories;  (1)  high  gradient  continua  -  including  higher  orders  of 
deformation  gradients,  and  (2)  first  gradient  continua  -  including  only  the  first  order  of 
deformation  gradients. 


Particle  spin  and  moment  transmitting  through  the  continua  are  considered.  The 
continuum  variables  involves  higher  orders  of  both  displacement  and  rotation  gradients. 

Class  2  :  High  Gradient  Couple  stress  /  Cosserat  Continua 

The  particle  spin  is  neglected,  thus  the  particle  rotation  is  induced  solely  from  the  rigid 
body  rotation  of  displacement  field.  The  gradients  of  rigid  body  rotation  can  still  cause  the 
transmitting  of  contact  moments  through  granular  material.  This  mechanism  resembles  that 
proposed  in  couple  stress  theory  and  Cosserat  theory  (Cosserat  1909,  Truesdell  and  Toupin 
1960,  Toupin  1962,  Mindlin  and  Tiersten  1962).  The  moment  transmitting  are  found  to  have 
interesting  effects  on  the  propagation  and  dispersion  of  waves  in  granular  media  (Chang  and 
Gao  1995b). 

Same  as  class  2,  we  further  neglect  the  moment  transmitting  caused  by  the  gradients 
of  rigid  body  rotation.  The  continuum  variables  still  involves  higher  orders  of  displacement 
gradients.  Chang  and  Gao  (1995a)  have  derived  this  class  of  constitutive  models  for  granular 
materials.' 
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Using  up  to  the  first-order  gradients  of  displacement  and  rotation.  The  continue 
resemble  the  micro-polar  type  proposed  by  Eringen  (1968),  Constitutive  law  of  this  type  for 
granular  material  shows  that  moment  trasmitting  can  effect  significantly  the  stress  distribution 
under  loads  (Chang  and  Ma  1991,  1992). 


Neglecting  the  particle  rotation  gradient,  the  particles  within  a  group  of  particles  spin 
with  the  same  magnitude.  In  this  situation,  the  granular  medium  does  not  transmit  moment 
thus  is  no  longer  a  micro-polar  medium.  However,  since  the  particle  spin  is  still  in  presence, 
the  granular  medium  is  of  a  quasi-micropolar  type.  Chang  and  Misra  (1989)  and  Chang 
(1993)  have  derived  this  class  of  constitutive  models  for  granular  materials. 

With  further  neglect  of  the  particle  spin  from  class  5,  the  particle  rotation  is  identical  to 
the  rigid  body  rotation.  Constitutive  law  for  granular  medium  thus  derived  resembles  that  of 
classic  continua  in  solid  mechanics.  Constitutive  models  of  this  class  for  granular  material 
can  be  found  in  the  work  by  Bathurst  and  Rothenberg  (1988),  Chang  (1988),  Jenkins  (1988), 
Walton  (1987),  etc.. 

Most  work  in  the  literature  on  micro-macro  properties  of  granular  material  treat  granular 
material  as  first  gradient  continua  in  the  classes  4-6.  Due  to  the  difficulties  of  dealing  with  the 
high-rank  tensors,  very  little  attention  has  been  paid  to  the  representation  of  granular  material 
as  high  gradient  continua  in  the  classes  1-3.  To  illustrate  the  relationship  between  micro¬ 
macro  properties  such  as  internal  length  of  granular  material  and  inter-particle  stiffness,  we 
discuss  the  constitutive  constants  of  a  high-gradient  model  for  isotropic  granular  materials. 

3.4  Constitutive  Relationship  for  an  assembly 


For  a  representative  volume  of  the  granular  material,  the  medium  can  be  treated  as 
statistically  homogeneous  and  can  thus  be  regarded  as  centrical  symmetry.  Under  this 
situation,  the  effect  of  the  first  gradient  of  strain  and  rotation  can  be  neglected  and  the 
constitutive  equation  for  the  granular  assembly  is  in  the  following  form 


This  type  of  high-gradient  constitutive  law  for  granular  material  is  in  the  same  class  of 
constitutive  law  proposed  by  Toupin  and  Gazis  (1964),  Mindlin  (1965),  Coleman  and 
Hodgdon  (1985),  etc. 

Consider  an  idealized  granular  material  which  is  made  of  equal-size  spheres  with  an 
isotropic  packing  structure;  The  material  has  identical  inter-particle  stiffness  for  all  contacts 
which  is  independent  of  contact  force.  The  expression  of  constitutive  relation  is  given  by 
(Chang  and  Gao  1995a) 


+C^  +C2(V®e  j 


(3.4) 


In  this  equation,  X  and  ji  are  the  Lame  constants  which  relate  to  the  inter-particle 
stiffness  by 


x=4a(/i:,-/g 

li=2a(2fi:,+3Jg 


(3.5) 


Mr 

where  represents  the  density  of  packing  structure.  M  is  the  total  number  of  inter¬ 


particle  contacts  in  the  representative  volume  V,  and  r  is  the  radius  of  particles.  The  other  two 
constants  corresponding  to  the  high-gradient  terms  are  as  follows: 


(3.6) 


3.4.1  Micro-macro  Properties 

Compare  the  derived  stress-strain  law  with  generalized  Hooke's  law,  the  constants  X 

and  (1  are  the  usual  Lame  constants.  The  corresponding  Young's  modulus  and  Poisson's 
ratio  are  derived  as: 

This  equation  provides  a  method  for  estimating  elastic  moduli  based  on  the  values  of  inter¬ 
particle  stiffness.  Eq.  3.5  can  also  be  rearranged,  given  by 
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(3.7) 
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(3.8) 


which  provi(des  a  means  for  estimating  inter-particle  stiffness  basecd  on  the  measured  value  of 
elastic  moduli. 

The  values  of  high  gradient  constitutive  coefficients  and  are  expressed  in 

terms  of  contact  stiffness  as  given  in  Eq.  3.6.  Substituting  Eq.  3.8  into  Eq.  3.6,  the  values  of 
and  C2  can  be  expressed  in  terms  of  Lame  constants  as  follows: 


3(A.-p)' 

10p  J 


or 


/ 
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(3.9) 


Eq.  3.9  provides  a  useful  method  for  estimating  the  high  gradient  constitutive  constants, 
.  directly  from  the  Lame  constants  and  particle  size.  This  relationship  is  of  practical 

use  since  the  high  gradient  constitutive  coefficients  are  difficult  to  be  evaluated  from 
laboratory  test. 


3.4.2  Role  of  Internal  Length  and  Inter-particle  Stiffness 

The  high-gradient  constitutive  constants,  C■^ ,  C2  are  functions  of  internal  length  r  of 

the  granular  assembly  which  represents  the  particle  size.  When  the  particle  size  of  a 
granular  material  is  very  small  compared  with  the  representative  volume  of  the  material,  the 
effect  of  strain  gradient  in  the  constitutive  equation  can  be  neglected.  Thus,  the  constitutive 
equation  reduces  to  the  generalized  Hooke’s  law  for  granular  materials.  On  the  other  hand, 
the  effect  of  strain  gradient  becomes  pronounced  as  the  particle  size  increases. 
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Inter-particle  stiffness  have  significant  effect  on  the  second-gradient  constitutive 
coefficients.  Three  cases  are  noted; 

(1)  kjk„  is  less  than  1;  This  corresponds  to  particles  with  smooth  surface.  Such  a 

ratio  of  inter-particle  stiffness  leads  to  a  positive  Poisson's  ratio  and  thus  a  positive  .  X  and 
a  positive  .  As  the  ratio  of  KjK„  decreases  (i.e.,  smoother  surface),  values  of  C, 

and  X  increase  while  value  of  Cj  decreases.  The  limiting  value  physically  possible  for 

KjK„  is  zero.  Corresponding  to  this  limiting  condition,  it  is  noted  that,  under  the  present 

formulations  for  granular  assemblies  with  spherical  particles,  the  predicted  Poisson’s  ratio  can 
not  be  greater  than  0.25  and  the  Lame  constant  X  can  not  be  greater  than  n  . 

(2)  KjK^  is  equal  to  1:  This  corresponds  to  particles  with  rough  surface.  Such  a 

ratio  of  inter-p4rticle  stiffness  leads  to  a  zero  Poisson’s  ratio  and  thus  a  zero  A,  and  a  zero 
.  Under  this  condition,  Cj  is  a  positive  number. 

(3)  KJK,  is  greater  than  1 ;  This  corresponds  to  particles  with  very  rough  surface, 

which  is  an  unusual  situation.  Such  a  ratio  of  inter-particle  stiffness  leads  to  a  negative 
Poisson’s  ratio  and  thus  a  negative  X  and  a  negative  . 

3.5  Comparison  with  Other  Models 

A  convenient  simple  form  based  on  the  second-gradient  theory  can  be  arrived  by 
neglecting  the  effects  of  volume  strain  gradient  in  Eq.  3.4,  and  by  assuming  the 

condition  of  Arj=0. 
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The  simplified  version  of  Eq.  3.4  takes  the  following  form  of  constitutive  equation 

(3.10) 


It  is  noted  that  the  simplified  model  of  second-  gradient  theory  (with  K^=  0)  given  in  Eq. 
3.10  is  similar  to  the  form  given  by  Bazant  and  Pijaudier-Cabot  (1989)  based  on  nonlocal 
continuum  theory,  and  the  linear  model  of  Mindlin's  second-gradient  theory  (Mindlin,  1965, 
Beran  and  McCoy,  1970b). 

Comparison  between  nonlocal  model  and  high-gradient  model  for  composite  media  can 
be  found  in  the  work  by  Beran  and  McCoy  (1970a,  b)  and  Levin  (1971),  Here  we  discuss  the 
relations  between  nonlocal  theory  and  the  present  high-gradient  model  for  granular  material. 
Based  on  Taylor’s  expansion  (Eq.  3.2),  the  high  gradients  of  displacement  at  the  local  point 
'n'  contribute  to  the  displacement  at  another  local  point 'm'  in  the  neighbouring  distance. 
Therefore,  second-order  gradient  terms  in  the  present  model  reflect  the  effect  of  neighbouring 
medium  on  a  local  point.  It  is  thus  expected  the  high-gradient  model  is  equivalent  to  the 
nonlocal  model  under  certain  conditions. 

In  the  Eringen-Kroener  model  of  nonlocal  elasticity,  the  constitutive  equation  is 


Off  W=/^a(l*  -Jf1)  5ff  ettCsO  ^2n  ( 

where  is  the  non-local  function  (see  Kroener,  1967;  Eringen  and  Edelen,  1972; 

Eringen,  1973).  Expand  the  strain  in  the  integrand  at  point  x,  i.e. 


and  substitute  the  series  of  the  strain  in  Eq.  3.12  into  Eq.  3.11  of  nonlocal  model.  After 
neglecting  higher  terms,  we  obtain  the  first-order  approximation  of  nonlocal  model,  given  by 


=  (1  +-cV“)  (A.  ji  tg) 

C  =f^^a{\x'-x\]  {xl,-xfdv(x) 


(3.13) 


in  which  the  integrand  of  Eq.  3.13  and  the  resultant  c  is  positive.  The  value  of  ^.depends  on 
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the  nonlocal  function.  Note  that  the  first-order  approximation  of  nonlocal  model  given  in  Eq. 
3.13  is  same  as  the  simplified  form  of  the  present  high-gradient  model  in  Eq.  3.10,  and  the 

p 

constant  in  Eq.  3.13,  c=^r^. 


3.6  Summary 

In  this  chapter,  classes  of  equivalent  continua  for  granular  materials  are  described.  As 
expected,  during  the  transformation  from  a  complicated  discrete  system  to  a  simpler 
continuum  system,  certain  information  are  lost.  Therefore,  the  constitutive  relations  for  the 
equivalent  continuum  can  not  respond,  on  the  macro  scale,  all  deformation  features  just  as 
the  original  granular  material  does.  How  simple  the  equivalent  continuum  should  be  depends 
on  how  much  “essential  behavior"  of  the  granular  materials  are  of  interest.  Unfortunately,  the 
limitations  of  modelling  capability  are  inherited  from  the  type  of  continuum.  Therefore  the 
required  modelling  capability  for  problems  encountered  determines  what  type  of  equivalent 
continuum  to  be  selected  for  analysis.  For  example,  the  continuum  of  micropolar  type  is 
useful  to  the  study  of  moment  transmitting  in  granular  media  while  the  high-gradient  continua 
is  useful  to  the  study  of  particle  size  effect  on  wave  dispersion. 
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CHAPTER  4 


Wave  Propagation  in  Granular  Material 
Using  High  Gradient  Theory 


Treating  granular  materials  as  a  high  gradient  continue  media,  the  propagation  of  wave 
in  granular  materials  are  analyzed.  The  constitutive  coefficients  of  the  high  gradient  continue 
are  expressed  in  explicit  functions  of  inter-particle  stiffness.  In  contrast  to  the  classical 
continue,  the  high  gradient  terms  represent  the  internal  characteristic  length  of  granular 
materials.  With  the  high  gradient  model,  we  solve  the  wave  propagation  in  a  finite  strip  of 
granular  material  using  the  method  of  separated  variables,  in  which  the  orthogonal 
eigenfunctions  can  be  obtained  from  the  boundary  conditions.  The  short  waves  with  lengths 
less  than  about  two  times  of  the  particle  diameter  can  not  propagate  through  the  granular 
material  because  of  the  discrete  nature  of  the  granular  material.  The  solution  derived  from  the 
high  gradient  theory  consists  of  a  series  of  harmonic  waves  which  are  admissible  of 
propagation  in  the  granular  material.  The  present  numerical  results  show  that  due  to  the 
nonlinear  dispersion  of  wave  in  granular  medium,  the  peak  value  of  the  stress  wave  will  be 
decay  as  wave  propagating.  A  series  smaller  stress  waves  including  extension  and 
compression  waves  will  be  emerge  after  the  main  stress  wave  goes  through,  which  is  caused 
by  the  reflection  of  the  stress  wave  on  boundary  of  inter-particles.  The  predicted  phenomena 
are  qualitatively  in  agreement  with  the  simulation  results  from  discrete  element  methods  and 
with  the  experimental  observations  from  a  chain  of  photo-elastic  disks  subjected  to  a  impact 
load. 


4.1  Introduction 

Dynamic  response  of  granular  material,  such  as  soil,  ceramic  and  powder  material,  has 
been  of  interest  to  the  field  of  civil  engineering,  material  processing  and  geophysics.  This 
chapter  is  aimed  to  study  wave  propagation  in  granular  material  treated  as  a  high  gradient 
media. 

The  high  gradient  elastic  stress-strain  model  used  in  this  chapter  was  developed  from  a 
micromechanics  approach,  in  which  the  constitutive  coefficients  are  expressed  as  explicit 
functions  of  inter-particle  stiffness  and  particle  size  (C.S.  Chang  and  J.  Gao,  1995a,b).  In  this 
chapter,  the  high  gradient  model  is  applied  to  analyze  the  wave  propagation  in  a  finite  strip 
fixed  at  one  end  and  subjected  to  dynamical  traction  at  the  other.  The  solution  cambe 
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obtained  from  the  solutions  of  a  series  of  harmonic  waves  using  the  method  of  separation  of 
variables.  The  solution  proves  that  the  short  waves  with  lengths  less  than  about  two  times  of 
the  particle  diameter  can  not  propagate  through  because  of  the  discrete  nature  of  the  granular 
material.  The  solution  also  shows  a  decay  of  the  peak  of  stress  wave  as  the  propagation  of  the 
wave.  This  phenomenon  is  similar  to  the  decay  of  peak  contact  force  observed  from 
experiments  by  A  Shukla  and  C.  Damania(1988)  and  from  Discrete  Element  Analysis  by  C. 
Thornton  and  Randall,  C.W.(1988),  Sadd  etc(1992),  J.S.  Lee  etc(1992),  for  a  chain  of  disks. 
This  chapter  discuss  the  decay  of  the  peak  stress  influenced  by  the  size  of  particles. 

In  this  microstructural  continuum  model,  the  stress  wave  excited  by  dynamic  impact 
consists  of  both  compression  and  extension  waves  propagating  along  the  granular  material. 
This  is  different  from  a  simple  compression  wave  in  classic  solution  of  elastic  material.  The 
phenomenon  of  tension  reflection  is  similar  to  that  from  Discrete  Element  analysis  (J.S.Lee  etc. 
1992).  The  solution  shows  that  this  high  gradient  continuum  model  captures  important  salient 
features  of  granular  material. 


4.2  Wave  Propagation  in  a  Finite  Strip  of  Granular  Materials 

Granular  material  can  be  treated  as  a  collection  of  particles.  The  geometric  deformation 
in  discrete  system  are  described  by  the  translations  «"  and  rotations  o"  (n=1,2,...).  To 
develop  a  continuum  mechanics  model  for  the  behavior  of  particle  assembly,  it  is  desirable  to 
view  the  discrete  system  as  an  equivalent  continuum  system.  To  this  end,  we  define  the 
macro-scale  continuum  geometric  deformation  fields  as:  displacement  Ufix)  and  rotation 
.  When  x  is  located  at  the  centroid  of  n-th  particle,  the  continuum  geometric  fields 

are  compatible  with  that  of  the  discrete  system,  i.e.  and  (iij[x’^  =  (a"  .Therefore 

the  inter-particle  displacement  and  inter-particle  rotation  can  be  expressed  as  a 

series  of  high  gradients  of  the  particle  deformations.  From  the  kinematic  hypothesis,  the 
constitutive  law  of  granular  material  with  random  packing  structure  can  be  expressed  as 
follows  (C.S.  Chang  and  J.  Gao  1995b): 

where  the  constitutive  constants  are  dependent  of  inter-particle  properties  and  the  size  of 
particles,  given  by 
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A=4a{X,-iy, 

c,=a,(/s:„-jg, 


,i=2«(2i;+3^,) 


C2=2a,{ji:,+|jg 

in  which  M  is  the  total  number  of  inter-particle  contacts  in  the  representative  volume  V;  r  is  the 
radius  of  particles;  and  [K^,  K^,  G^,  GJ  are  the  distortion  and  rotation  contact  stiffness 
constants  along  normal  and  tangential  directions  respectively. 

We  now  consider  the  case  of  wave  propagation  in  a  finite  strip  of  granular  material.  In 
this  case,  only  the  displacement  i£(jc,i)  and  stress  in  x-direction  are  considered. 

The  constitutive  equation  of  stress  o(*,^  in  the  aforementioned  high  gradient  model  is  given 
by 

(4.2) 

dx 


where 


X42^  7 


With  the  high  gradient  model,  the  one  dimensional  wave  equation  becomes 


P  ^  =i£ =(;t42n)(1  +c-^)-^ 
*2  ar  ^  ac2  3*2 


(4.3) 


where  p  is  the  material  mass  density. 

It  is  noted  that  the  wave  equation  is  a  fourth-order  differential  equation.  When  the 
influence  of  the  internal  length  on  material  behavior  is  neglected,  i.e.  c=0,  the  wave  equation 
in  Eq.  4.3  reduces  to  the  classical  form  in  elastic-dynamics.  The  wave-propagation  in  the 
finite  strip  requires  four  boundary  conditions  and  two  initial  conditions.  Two  of  the  boundary 
conditions  are  given  by 


«(0,#)=0:  =  l>0  (4.4) 

3*  >.+2|i 

The  two  extra  boundary  conditions,  by  use  of  variational  principle,  should  be 

Mh91^Q  (4.5) 

3*2  ’  3*2 


f)  -  Q.  3^Bfe  0)  _  Q 
3*2  ’  3t? 
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The  initial  conditions  are 


k(x:,0)  =  Wo(x);  =  (KpSL  (4.6) 

at 

In  the  following,  two  types  of  excitation  are  discussed,  namely,  initial  disturbance  and 
dynamic  traction. 


4.3  Wave  Propagation  Subject  to  Initial  Disturbances 


First,  we  consider  that  the  wave  in  the  finite  strip  is  excited  by  the  initial  disturbance 
which  described  by  the  initial  conditions  Eq.  4.6  with  the  homogeneous  boundary  conditions. 

To  obtain  the  solution  by  the  method  of  eigenfunction  expansion,  we  further  propose 
the  solution  to  be  in  the  form  .  Substituting  the  assumed  solution  into  Eq.  4.3 

and  separating  the  variables,  we  can  express  the  wave  equation  as  two  ordinary  differential 
equations,  given  by 


cX^(*:)4-X"{x)+-^Jir=0 


(4.7) 


with  the  boundary  conditions 


x(0)=o,  x"(0)=0:  r(L)=o,  x"%L)=Q 


(4.8) 


Clearly,  the  function  8Xp(ic»i)  is  a  solution  of  Eq.  4.7  which  corresponds  to  the 
harmonic  wave  propagating  in  material.  The  constant  <ii  is  the  natural  frequency  of  the 
system.  The  solution  of  equation  (4.7b)  leads  to  a  series  of  eigenvalues  (n  =  1,2,3,...)  .  A 

particular  mode  of  the  oscillation  is  termed  the  eigenfunction  associated  with  the 


f 

natural  frequency  .  The  eigenvalues  satisfy  the  following  theorems: 


Ci3„ 

Theorem  1.  All  eigenvalues  P,=—  of  Eq.  4.7  are  real  when  the  natural  boundary 

Y2 

conditions  are  satisfied. 


Theorem  2.  The  eigenfunctions  associated  with  real  p,  constitute  an  orthogonal  set.  • 
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It  is  noted  that  the  solution  of  Eq.  4.7b  is  dependent  on  the  constitutive  parameter  of 
high  gradient  term. 


(a)  When  4cp>  1  ,  the  solution  can  be  expressed  as  follows 


X{x)  =ex.p{-K^x)\ACQS{K2x)  +  JJsinCiTjJr)]  +exp(i:^j:)[Ccos(A:a;c)  +jDsin{Ji:;jai:)]  (4.9) 


where  K. 


''-i 


c  2 


•&oos|. 

c  2 


To  satisfy  the  boundary  conditions  given  by  Eq.  4.8,  we  substitute  Eq.  4.9  into  Eq.  4.8 
leading  to  a  set  of  linear  algebra  equations  with  the  coefficients  A,B,C  and  D.  It  can  be  seen 
that  non-trivial  solution  of  A,  B,  C  and  D  can  be  obtained  if  and  only  if  the  parameter  matrix  is 
zero.  Since  the  parameter  matrix  of  the  linear  equations  is  not  equal  to  zero,  i.e. 


2K^K2{Ki  +  x|)(cos^(JSy:)  ch^(K^L) + sm^{K2L)sh^{K^L))  *Q  (4  1 0) 


Therefore,  the  constants  A  =  B  =  C  =  D  =  0. 

(b)  When  4cp  s  1  ,  the  solution  of  eigenfunction  can  be  expressed  as  the  following 

form 


X{pc)  =A  cos(— +v/1  -4cp)  +  B  sin(— +v/l  -4<;p) 

VS  ^ 

*  ccos(— *  o  sin(— 

V®  V2c 


(4.11) 


From  the  boundary  conditions,  we  can  again  obtain  the  linear  algebra  equation  of  the 
constants  A,  B,C  and  D.  At  jc=0  ,  the  linear  algebra  equations  gives  A=C=0  .Applying 
the  boundary  conditions  on  another  end  of  the  strip,  it  can  be  seen  that  non-trivial  solution  of  B 
and  D  can  be  obtained  if  and  only  if  the  parameter  matrix  is  zero,  i.e. 


Cos(-^^/nVT^)  =  0  (4.1 2) 

VS 


cos(— /rVT^)  =  0 


In  fact,  the  two  cases  lead  to  the  same  conclusion  which  gives 
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= in  *■  hr,  (« = 1 .2,3...) 

From  Eq.  4.1 1  we  can  obtain  the  eigenvalues  given  by 

P, - ^(1  -  (n=1,2.....) 


(4.13) 


(4.14) 


Because  the  restriction  condition  0<4cp,<  1  ,  the  eigenvalues  can  only  be  a  finite 
series  instead  of  an  infinite  series.  The  number  of  terms  of  the  series  is  determined  from  the 
relation  . 


N  =  max. 


(4.15) 


A  finite  series  of  proper  orthogonal  eigenfunctions  exist  in  the  high  gradient  theory, 
given  by 

1 


XJliK)  =(|)2sin[(«^l)|x|;  «  =  0,1,2.....]y. 


(4.16) 


From  the  above  discussion,  the  waves  smaller  than  certain  wavelength  can  not  exist  in 
the  high  gradient  medium  due  to  the  effect  of  internal  length.  The  same  conclusion  was 
derived  from  the  dispersion  relation  of  granular  media  (Chang  and  Gao  1995b). 

The  max'imum  wave  number  and  the  minimum  wavelength  can  be  obtained  from  the 
restriction  of  non-negative  eigenvalues.  From  Eq.  4.12,  we  obtain  the  maximum  wave  number 


or  the  minimum  wavelength 


^mln  "  .  “  2nVc  •=  Cff 

*mex 


(4.17) 


(4.18) 


where  =  3.35 


N 


is  a  constant. 


For  the  limiting  case  of  ir^=0  ,  the  minimum  wavelength  is  3.45r,  approximately  two 
times  of  particle  diameter,  below  which  the  waves  can  not  propagate  in  the  finite  strip. 

From  the  wave  numbers  obtained,  we  can  deduce  the  natural  frequencies  given  by 


1_4(„4)V; 


»  =  1i2 . N. 


(4.19) 


L  N  2! 

Corresponding  to  natural  frequencies,  the  solution  of  Eq.  4.7a  can  be  in  thg  form 
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n=0,i,2,...,N. 


(4.20) 


r,=a,cos{(o,0  +fe„sin(w,/); 


Thus  the  solution  of  the  wave  propagation  in  the  finite  strip  with  a  fixed  end  and  free 
end  is  expressed  as  follows 

«(r,  /)=£  {a,C0S(£J^)  +  &,Sin(co^))sinI^{n+-l)n]  (4  2 

where  coefficients  a^,  are  to  be  determined  from  the  initial  conditions,  given  by 


|/Q\Ws'n[y(n+^)iiI  dX‘, 


(4.22) 


4.4  Wave  Propagation  Subject  to  Dynamic  Traction 

In  this  section,  we  will  consider  wave  propagation  in  the  finite  strip  excited  by  the 
dynamic  traction  given  by  the  boundary  conditions  Eqs.  4.4  and  4.5,  with  zero  initial 
conditions.  In  this  case,  we  propose  the  solution  to  be  in  the  form 


k(x,i)=}^,i) 

0 


(4.23) 


where  are  the  eigenfunctions  given  by  Eq.  4.16;  and 

are  the  unknown  functions. 

Substituting  Eq.  4.23  into  the  wave  equation  (4.3)  leads  to  the  following  two  differential 
equations; 


Q  dr 


_Q 

at*  at* 


(4.24) 


Next,  multiplying  both  sides  of  Eq.  4.24a  by  XJ^)  ,  integrating  over  [0,L],  and  from 
the  orthogonal  property  of  eigenfunctions,  we  have 

^{f)4.«X(/)  =  4.^(f)  (4  25) 

where 


<b„{t)  =  -f^¥ipc,t)X,(^)dx 
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For  zero  initial  conditions,  we  obtain  the  solution  of  Eq.  4.25,  given  by 


Sin((D,(/-T))dT  (4.26) 

It  is  noted  from  Eq.  4.25b  that  we  can  approximately  express  the  function 

0 

Substituting  Eq.4.26  into  Eq.  4,23  leads  to  the  following  form  of  the  solution 

=  "E  sin{a,{/-T))dT  (4.27) 

In  addition,  the  unknown  function  is  determined  from  the  static  equation 

(4.24b).  Its  solution  is 

Y{x,t)  =  a(t)  +  b{t)x  +  c(f)sin(— )  4  d{t)  COS(— )  (4.28) 

}/c  ^ic 

From  the  boundary  conditions  (4a), (5a)  and  (6b),  we  obtain  that  a=c=d=0.  And  then, 
from  Eq.  4.4b,  the  solution  of  Eq.  4.24  becomes 


i^,r)  = 


k  +2n 


(4.29) 


Substituting  Eq.  4.29  into  Eq.  4.25b,  we  obtain 

(4.3 

And  then,  substituting  Eq.  4.30  into  Eq.  4.27  leads  to  the  final  solution  of  wave  propagation 
subject  to  dynamic  traction,  given  by 


{X  +2ii)n^  Q  {2n+1)''  L  Z  •>» 


For  a  specific  time-dependent  traction  p{H  given  by 

p(0=Po: 

=  0;  tQ<t 


(4.32) 
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The  solution  can  be  obtained  by  substituting  Eq.  4.32  into  Eq.  4.31,  thus 


(A.  +2(1)31  Y  (2n+1)2  L  2 

And  then,  the  stress  field  is  derived  from  the  constitutive  equation  (4.2),  given  by 


(4.33) 


t)  =  8PoX:  ~(«^|)"n')CX)S(^(«+|)lT)/« 


(4.34) 


where 


/W=^[1-cos(«,r)): 

=  sin((o„|)sin(«,(#-|)); 


0<f<<o 


4.5  Summary 

From  Eq.  4.18,  the  admissible  minimum  wave  length  is  proportional  to  the  particle  size. 
The  number  of  terms  of  the  harmonic  waves  propagating  in  granular  materials  is  inversely 
proportional  to -the  particle  size.  When  the  diameter  of  particles  is  very  small  compared  to  the  , 
the  effects  of  high  gradient  terms  in  the  constitutive  equation  are  negligible,  and  the  solution  of 
wave  propagation  is  degenerated  into  classical  solution. 

In  classical  solution,  all  harmonic  waves  have  the  same  wave  velocity.  In  high  gradient 
theory,  each  of  harmonic  wave  has  a  different  velocity.  Since  the  excited  wave  due  to  dynamic 
traction  is  a  summation  of  harmonic  waves,  the  shape  of  the  excited  wave  keeps  constant 
during  propagating  in  the  classical  theory.  However,  in  high  gradient  theory,  the  shape  of  the 
excited  wave  changes  due  to  dispersion. 

Eq.  4.34  is  used  for  analysing  the  stress  wave  propagating  along  the  strip.  A  typical 
sequence  of  the  stress  wave  propagation  in  the  finite  strip  subjected  to  dynamic  traction  is 
shown  in  Fig.  4.1  using  dimensionless  variables. 

Although  the  stress  wave  excited  at  the  beginning  of  dynamic  impact  is  only  a  simple 
compression  wave,  at  later  stage,  the  stress  waves  include  both  compression  waves  and 
extension  waves.  The  phenomena  is  caused  by  the  dispersion  of  stress  waves  due  to  the 
discrete  nature  of  granular  material. 

The  normalized  peak  of  stress  wave  decay  with  time.  The  rate  of  decay  decrease  with 
time.  The  influence  of  particle  size  on  the  decay  stress  wave  during  propagation  is  shown  in 
Fig.  4.2.  When  the  particle  diameter  increase,  the  phenomenon  of  decay  is  a  significant  factor 
in  the  analysis  of  wave  propagation  in  granular  media.  - 
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CHAPTERS 

Normal  and  Tangential  Compijance  for 
Elastic  Conforming  Binder  Contact. 


This  chapter  describes  the  normal  and  tangential  conforming  contact  compliance  for  a 
system  of  two  elastic  particles  bonded  by  a  layer  of  elastic  binder  in  between.  The  governing 
equation  of  this  problem  is  a  second  kind  of  Frediholm  integral  equation  with  singularities  of 
logarithmic  type.  Exact  solution  for  the  unknown  interfacial  pressure  between  particle  and 
binder  is  difficult  to  arrive.  Derivations  of  compliance  are  presented  in  the  forms  of  the  upper 
and  lower  bounds,  and  of  the  best  estimate  based  on  physical  approximations.  It  shows  that 
the  derived  elastic  compliances  agree  favorably  with  those  of  "discretized  exact  solutions" 
obtained  from  numerical  methods. 

5.1  Introduction 

The  subject  of  layer/binder  contact  frequently  occurs  in  granular/particulate  materials 
such  as  asphaltic  concrete  or  cemented  sand.  This  subject  is  also  important  in  tribology, 
involving  the  mechanical  behavior  of  coated  materials.  Many  topics  in  this  area  have  been 
investigated  in  the  past  years.  For  example,  Muki  (1960)  studied  the  problem  of  contact 
between  a  layer  and  an  elastic  half  space.  Goodman  and  Keer  (1975)  studied  a  case  of 
surface  layers  bounded  to  a  substrate.  Bentall  and  Johnson  (1968)  worked  on  a  plane  strain 
layered  problem  which  was  further  studied  by  Meijers  (1968)  and  Alblas  and  Kuipers  (1970) 
for  both  conditions  of  thin  and  thick  layers.  Matthewson  (1981)  presented  a  theory  of 
indentation  of  a  soft  thin  coating  by  a  rigid  body.  Keer  et  al.  (1991)  investigated  the 
compliance  of  coated  elastic  bodies  in  contact.  Dvorkin  et  al.  (1991)  employed  numerical 
solutions  to  examine  the  normal  interaction  of  two  elastic  spheres  separated  by  an  elastic 
cementation  layer  and  recently  (1994)  extended  the  numerical  solutions  to  examine  tangential 
deformation  of  two  cemented  spheres.  In  addition,  many  related  topics  can  be  found  in  the 
books  by  Johnson  (1985)  and  by  Gladwell  (1980). 

This  study  is  focused  on  the  compliance  of  a  system  of  two  elastic  particles  bonded 
together  by  a  thin  layer  of  binder.  We  aim  to  derive  closed-form  relationships  between  the 
forces  and  the  relative  particle/binder  movements  in  this  system.  Closed-form  expressions  are 
of  particular  interest  because  they  can  be  readily  incorporated  into  discrete  element  methods 
for  the  analysis  of  a  large  assembly  of  particles. 

Progression  of  the  chapter  begins  with  establishing  integral  equations  that  govern  the 
interfacial  contact  pressure  distribution  between  particle  and  binder.  Analyses  of  the  upper  _ 
and  lower  bounds  with  respect  to  the  compliance  of  the  two  particle  system  are  followed.  The 
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best  estimate  of  compliance  based  on  physical  approximation  is  then  pursued.  The  closed- 
form  analytical  expressions  are  compared  with  the  numerical  solutions  obtained  directly  from 
solving  the  governing  integral  equations  by  a  discretization  technique. 

5.2  Formulation  of  the  Problem 


Fig.  5.1  shows  an  axi-symmetric  configuration  of  two  particles  bonded  by  a  binder  in  a 
cylindrical  coordinate.  The  function  z  =  h(r)  represents  the  geometry  of  interfacial  boundary 
between  the  particles  and  the  binder,  given  by 


h{r)=h^ 


1  — - 


(5.1) 


where  a  is  the  radius  of  contact  area,  ho  is  the  thickness  of  the  binder  at  r  =  0,  and  d  is  the 
dimensionless  shape  parameter  related  to  the  curvature  of  particle  surface,  which  is  limited  in 
a  range  0  <  d  <  1 .  For  a  planer  surface,  d  is  zero.  For  a  spherical  particle,  d  is  given  by 


d= 


(5.2) 


where  Fi  is  the  radius  of  the  spherical  particles. 

We  denote  the  constraint  modulus  E,  and  E2,  Poisson’s  ratio  u,  and  for  the  particles 
and  the  binder  respectively,  where  the  constraint  modulus  E,  and  Ej  are  defined  as 


'  1-2vj  ’ 


i=1,2 


(5.3) 


and  G,  and  Gj  are  the  shear  modulus  of  the  particles  and  the  binder  respectively. 

We  intend  to  derive  the  normal  and  tangential  compliance  of  this  two  particle  system 
with  an  elastic  binder.  The  governing  equations  that  relate  force  and  relative  movement  of  two 
particles  are  discussed  separately  for  the  normal  mode  and  the  tangential  mode. 


5.2.1  Normal  Compliance 

The  relative  normal  approach  5^  for  the  two  contact  bodies  is  separated  into  two 
components:  the  normal  displacement  at  the  binder-particle  interface  relative  to  the  particle’s 
centroid,  Wi(r);  and  the  normal  displacement  at  the  binder-particle  interface  (i.e.,  at  z=h(r)) 
relative  to  the  z=0  plane,  W2(r),  given  by 

Since  z  =  0  is  a  plane  of  symmetry,  the  binder  normal  displacement  vanishes  at  z=0.  As  the 
binder  is  a  thin  layer,  we  approximate  the  normal  strain  to  be  uniform  in  the  z  direction  across 
the  binder.  Thus  the  normal  displacement  W2(r)  contributed  from  binder  can  be  expressed  'as 
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follows: 


iV2(r)=ft(r)^  (5.5) 

£■2 

where  p(r)  is  the  interfacial  normal  pressure  between  the  particles  and  the  binder. 

In  the  analysis  of  the  normal  displacement  Wi(r),  we  assume  that  the  characteristic 
dimension  of  the  particle  is  much  larger  than  that  of  the  particle-binder  contact  area.  Thus  it  is 
justifiable  to  pursue  the  analysis  of  Wi(r)  based  on  a  half-space  premise.  Following  the  well- 
known  Boussinesq’s  equation,  Wi(r)  can  be  related  to  p(r)  by: 

2 

14,  (r)  =  ~  p(p)  P^P  (5  6) 

°  ®  r  -  2r  p  COS  0 

Substituting  Eqs.  5.5  and  5.6  into  Eq.  5.4,  we  have 


9  .  f‘ 

^  nBl  Jo 


where  l(p,  r)  are  defined  as 


7(p,r)=:/(t)  =  ;  k  =  (5 

V1  -  kCOSQ  r  +  P 

Integration  of  the  interfacial  pressure  function,  p(r),  over  the  contact  area  gives  the 


2  2 
n  +  p 


resultant  normal  contact  force  P, 


P,  =  2ix  r  p(r)r  dr 


Eqs.  5.7  and  5.9  indirectly  provide  the  compliance  relationship  between  the  relative  normal 
approach  ^  and  the  contact  force  through  the  interfacial  pressure  function. 


5.2.2  Tangential  Compliance 

The  relative  tangential  approach  in  the  x-direction  for  the  two  contact  bodies  is  also 
separated  into  two  components:  the  tangential  displacement  at  the  binder-particle  interface 
relative  to  the  particle’s  centroid,  u,(r,  0);  and  the  tangential  displacement  at  the  binder- 
particle  interface  (i.e.,  at  z  =  h(r))  relative  to  the  z  =  0  plane,  UgCr,  0),  given  by 

6,  -  Bi(r,e)  +  «2(r,e)  (5.’0) 

The  tangential  displacement  vanishes  at  the  plane  of  symmetry  z  =  0.  Considering 
that  the  binder  is  a  thin  layer,  we  approximate  the  tangential  strain  to  be  uniform  in  the  z 
direction  across  the  binder,  and  the  following  relation  can  then  be  derived: 

(5.11) 

G2 

where  G2  is  the  shear  modulus  of  binder,  q(r,  0)  is  the  interfacial  tangential  pressure  between 
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the  particles  and  the  binder. 

-  Since  the  characteristic  dimension  for  the  particle  is  much  larger  than  that  of  the 
contact  area,  we  use  the  known  relationship  between  u,(r,  G)  and  q(r,  0)  based  on  the  half¬ 
space  premise  (Johnson  1985).  Thus  the  tangential  displacement  for  the  two  contact 
bodies  is  now  equal  to  the  summation  of  Ui(r,  0)  and  U2(r,  0)  with  an  error  of  the  order  (u,)^ 
(Dvorkin  et  al.,  1994) 

6^  =  Hr,p,Q,^,v^)pd4>dp  (5.12) 

where 


1  "Vj 

F(r,p,0,(l),Vi)  =  — 


+  Vh 


^2 = (rcose  -  p  cos^f + (rsin0  -  p  sin(l))^ 


(ytX)S0-pCQS<l))2 


} 


(5.13) 


where,  and  are  respectively  the  shear  modulus  and  Poisson’s  ratio  of  the  particle. 
Integration  of  the  interfacial  pressure  function,  q(r,  0),  over  the  contact  area  gives  the  resultant 
tangential  contact  force 

^^9  dr  (5.14) 

Again  the  governing  equations  (12)  and  (14)  provide  the  compliance  relationship  between  the 
relative  tangential  approach  and  the  contact  force  P,,  through  the  unknown  interfacial 
pressure  function,  q(r,  0). 

In  fact,  the  interfacial  pressure  functions,  q(r,  0)  or  p(r),  can  be  determined  by 
simultaneously  solving  Eqs.  5.12  and  5.14  or  Eqs.  5.7  and  5.9.  Unfortunately,  the  governing 
equations  (7)  and  (12)  are  the  second  kind  of  Fredholm  integral  equations  with  kernels  which 
have  logarithmic  singularities.  To  this  type  of  integral  equations,  analytical  solutions  are  difficult 
to  arrive.  However,  with  special  care  on  the  singularities,  they  can  be  solved  using  a  numerical 
discretization  technique.  Details  can  be  found  in  Zhu  et  al  (1995). 


5.3  Solutions  for  Two  Extreme  Gases 

The  analytical  solutions  of  the  interfacial  pressures  p(r)  in  Eq.  5.7  and  q(r,  0 )  in  Eq. 

5.12  are  known  for  two  extreme  cases,  namely,  (1)  rigid  particle  case  (i.e.,  E,  -^ooand  G,  ->oo 
while  Ej  and  G2  are  finite),  and  (2)  rigid  binder  case  (i.e.,  E^  and  Gj  are  finite  while  Eg  — >00  and 
Gi  -^00).  The  compliance  relationship  under  these  two  extreme  conditions  are  described  in 
this  section. 
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Rigid  Particle  Case 

In  the  rigid  particle  case,  the  relative  movement  of  the  two  contact  bodies  is 
contributed  only  from  the  deformation  of  binder.  Thus 


6,  =  A(r)  ^ 


6,  =  ftW 


<7(r,e) 

Go 


(5.15) 

(5.16) 


Subsequently,  for  the  rigid  particle  case,  the  corresponding  normal  interfacial  pressure 
denoted  as  p,(r)  is  given  by 

Pz  Ih 


Pi(r)  '■ 


ica' 


Hr)X 


(5.17) 


(5.18) 


and  the  interfacial  tangential  pressure  q(r,  0)  becomes  independent  of  the  variable  0 
(denoted  as  qXO)  arid  reads 

^^(r)  = 

na^A(r)X 

where 

^  =  InQ  (5.19) 

d 

and  d  is  the  shape  parameter  defined  in  Eq.  5.1. 

Thus  the  normal  and  tangential  compliance  relationships  between  the  contact  force 
and  the  relative  approach  5^  become 

*0 


~  ’  ^^x~ 


na^E^X 

Tza^GoX 


(5.20) 


(5.21) 


Rigid  Binder  Case 

In  the  rigid  binder  case,  which  represents  the  well  known  rigid  punch  problem,  the 
normal  interfacial  pressure,  denoted  as  P2(r),  is  given  by 

(^22) 

and  the  tangential  pressure,  q(r,  0)  is  again  dependent  only  on  r  (denoted  as  q2(r)): 


^2^  = 


2ua 


(5.23) 
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For  this  case,  the  normal  and  tangential  compliances  are: 


; 

"  ^Zx^x  » 


2aE.^ 

(5.24) 

2-v, 

(5.25) 

It  can  be  easily  verified  that  under  both  extreme  conditions,  r  Pi(r),  r  PjCr),  r  q/r),  and 
rq2(r),  are  monotonically  decreasing  functions  in  the  range  of  0  <  r  <  a.  Indeed,  the  functions 
r  p(r)  and  r  q(r)  are  monotonic  for  any  pairs  of  E,  ,  Eg  and  G,  ,  Gj  as  verified  from  the 
numerical  method  given  in  Appendix-A. 


5.4  Upper  Bound  Solution 


Explicit  compliance  relationships  are  easily  derived  for  the  two  extreme  conditions. 
However,  for  general  conditions,  analytical  solutions  to  Eq.  5.7  and  Eq.  5.12  are  difficult  to 
obtain.  The  upper  and  lower  bounds,  and  the  best  estimated  solution  based  on  the  physical 
approximations  have  been  derived  by  Zhu  et.  al  (1995).  A  brief  summary  is  given  in  this 
section. 

The  governing  equations  are  first  transformed  to  be  suitable  for  the  use  of 
Chebyshev's  inequality  for  integral.  Using  the  principles  of  Chebyshev’s  inequality  for  integral, 
the  equations  can  be  greatly  simplified  thus  leading  to  the  upper  bound  solution  for  the 
normal  compliance 

a,  i  (5.26) 

where 

(5.27) 

ff(d)=5.699-2.404<f+1.495rf^-1.079j3+0.841  j4-0.689^  (5  28) 

1  +d 


and  the  upper  bound  solution  for  the  tangential  compliance: 


2bi 

5,  ^  (Cu  * 


where,  b^  is  defined  in  Eq.  5.27. 


(5.29) 
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5.5  Lower  Bound  Solution 

Similarly,  application  of  the  Chebyshev's  inequality  for  integral  results  in  the  lower 
bound  solution  for  the  normal  compliance: 

6,  k  (J-2C,,+C2,)P,  (6  30) 

where 

d2=-(1+0.5d)X<  1  (5.31) 

4 

and  the  lower  bound  solution  for  the  tangential  compliance: 

6  (6.32) 

Based  on  the  upper  and  lower  bound  solutions,  the  true  normal  compliance  must  be 
between  (0,2  +  b^  C2z)  and  (ba  +  Cjz),  and  the  true  tangential  compliance  must  be  between 
(C,x  +  b,  Cjx)  and  (bj  +  ^2>d-  Examining  the  range  of  values  of  b,  and  bj,  for  versus  the 
shape  parameters,  the  maximum  relative  difference  between  the  upper  and  lower  bounds  of 
compliance  is  from  18-20  %. 

5.6  Best  Estimate  Based  on  Physical  approximations 

Based*  on  physical  approximation,  it  leads  to  the  following  best  estimate  for  normal 
and  tangential  compliance  relationship: 

Eqs.  5.33  and  5.34  satisfy  the  two  extreme  cases:  (1),  rigid  particle  case  (E,  — >ooand  Eg  finite); 
and  (2),  rigid  binder  case  (E,  finite  and  Eg  — >(X)).  In  addition,  the  best  estimated  compliances 
fall  in  between  the  upper  and  lower  bounds,  i.e.,  the  following  inequalities  are  always  satisfied: 

62C1 +C2  <  q  +c^<  q  +b^  C2  (5  35) 

hCu*C^  <  (6-36) 

^  Vi 

Eqs.  5.33  and  5.34  indicate  that  the  overall  compliance  corresponds  to  a  serial 
connection  of  the  two  compliances  C,  and  Cg ,  where  C,  represents  the  compliance  of  particle 
and  Cg  represents  the  compliance  of  binder.  The  best  estimated  analytical  compliances  are  in 
good  agreement  with  the  compliances  numerically  calculated  (Zhu  et  al  1995), 


(5.33) 

(5.34) 
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5.7  Summary 

For  the  compliance  of  an  elastic  particle-binder  system,  the  governing  equation  is  a 
second  kind  of  Fredholm  integral  equation  with  singularities  of  logarithmic  type.  Exact  solution 
is  difficult  to  arrive.  Even  solved  numerically,  special  care  needs  to  be  taken  for  the 
convergency  problems  associated  with  the  singularities.  This  approach  makes  a  use  of  the 
principles  of  Chebyshev's  inequality  for  integral  which  allows  simplifications  of  the  governing 
equations  thus  yields  remarkably  simple  closed-form  expressions  for  the  upper  and  lower 
bound  solutions  of  com.pliances.  The  best  estimates  of  compliances  have  been  found  to  be 
favorably  in  agreement  with  the  numerical  solutions  and  thus  can  be  considered  as  a  good 
approximate  solution  for  this  type  of  problem. 
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CHAPTER  6 

SUMMARY  AND  CONCLUSION 


This  report  reviews  the  results  of  this  project  in  four  different  areas:  1)  descriptions  of 
micro-structure,  2)  micro-macro  relationship,  3)  classes  of  micromechanics  constitutive  theory, 
and  4)  contact  law  of  the  inter-particle  binder.  These  four  areas  are  essential  to  the' 
construction  of  a  micromechanics  theory  for  particulate  media. 

Chapter  two  describes  several  methods  of  micro-structure  characterization  and  how 
the  microstructure  description  is  integrated  in  the  formulation  to  relate  the  micro-macro 
mechanical  behavior.  The  discrete  variables  selected  in  the  micromechanics  model  depend 
greatly  on  the  type  of  micro-structure  descriptor.  Thus  the  micro-structure  descriptor  is  the  key 
element  in  micromechanics  theory  for  granular  material.  The  cell  model,  as  oppose  to  the 
branch  model,  is  found  to  be  useful  in  describing  spacial  variation  of  microstructure  (i.e.,  the 
heterogeneity  of  material).  However,  there  are  still  many  fundamental  issues  need  to  be 
resolved  on  the  mechanics  of  cell  interaction. 

Chapter  three  addresses  the  limitation  of  classic  continuum  in  modelling  granular 
media.  In  order  to  model  some  features  due  to  the  effect  of  microstructure  and  discrete  nature 
of  granular  particles,  generalized  continue  different  from  the  classic  continuum  are  necessary 
to  be  considered.  Several  new  classes  of  generalized  continue,  ranging  from  simple  to 
complex,  are  derived  for  representing  granular  media.  Degree  of  complexity  of  the  continuum 
is  proportional  to  the  modelling  capability.  For  example  the  micro-polar  continua  is  useful  in 
modelling  the  effect  of  particle  rotation.  The  high  gradient  model  is  useful  in  modelling  large 
deformation  especially  after  peak  strength  of  the  granular  material  and  in  modelling  some 
distinct  features  for  dynamic  wave  propagating  through  granular  media. 

In  Chapter  four,  the  high  gradient  model  is  used  as  an  example  to  analyze  the  wave 
dispersion  in  granular  media,  a  phenomenon  that  can  not  be  modelled  by  the  classic 
continua.  The  results  show  that  waves  with  different  frequencies  have  different  velocities  in  a 
granular  medium.  Short  waves  propagate  slower.  Very  short  waves  can  not  pass  through  the 
medium.  Subsequently,  a  relationship  is  derived  for  the  decay  of  peak  stress  wave  as  a 
function  of  particle  size.  This  phenomenon  is  very  different  from  the  theory  of  wave 
propagating  in  classic  continuum. 

Chapter  five  presents  part  of  the  results  from  the  collaborated  work  with  Wright 
laboratory  at  Tyndall  Air  Force  Base.  This  chapter  describes  the  normal  and  tangential 
conforming  contact  compliance  for  a  system  of  two  elastic  particles  bonded  by  a  layer  of 
elastic  binder  in  between.  The  principle  of  Chebyshev’s  inequality  is  found  to  be  a  useful 
concept  for  obtaining  solutions  for  the  governing  equations  of  Frediholm  type.  Continuation  of 
the  work  are  undertaken  for  studying  compliance  of  a  binder/particle  system  in  rolling  and 
twisting  modes. 
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